{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "#### Ornstein-Uhlenbeck process without parameter estimation\n", "\n", "The Ornstein-Uhlenbeck (OU) process is a continuous-time continuous-space\n", "stochastic process governed by the stochastic differential equation (SDE)\n", "\n", " dY(t) = (Θ_1 - Θ_2 Y(t)) dt + Θ_3 dW(t)\n", "\n", "where W(t) is the Wiener process (standard Brownian motion). This OU process is\n", "stationary, Gaussian and Markovian. It is mean-reverting in the sense that it\n", "drifts towards its long-term mean of Θ_1 / Θ_2.\n", "\n", "For a step size Δt, the transition distribution can be obtained analytically as\n", "a Normal with mean μ(Δt) and variance σ^2(Δt):\n", "\n", " Y(t + Δt) | Y(t) ~ Normal(μ(Δt), σ^2(Δt))\n", " μ(Δt) = Θ_1 / Θ_2 + (Y(t) - Θ_1 / Θ_2) exp(- Θ_2 Δt)\n", " σ^2(Δt) = Θ_3^2 / (2 * Θ_2) * (1 - exp(- 2 Θ_2 Δt))\n", "\n", "Given a sequence of noisy observations of the process Y(t_n) ≈ y_n for\n", "1 <= n <= N, the task is to infer the posterior distribution over each Y(t_n).\n", "This choice is for the sake of simplicity, the example can be modified to\n", "infer posterior distributions over any Y(t).\n", "\n", "This example has been inspired by the LibBi example at\n", "http://libbi.org/packages/OrnsteinUhlenbeckBridge.html.\n", "\n", "##### References:\n", "- Y. Aït-Sahalia. Transition densities for interest rate and other nonlinear\n", " diffusions. The Journal of Finance, 54(4):1361–1395, 1999. ISSN 1540-6261.\n", " doi: 10.1111/0022-1082.00149.\n", "- L. Sun, C. Lee, and J. A. Hoeting. Penalized importance sampling for parameter\n", " estimation in stochastic differential equations. 2013.\n", " URL http://arxiv.org/abs/1305.4390.\n", "- Del Moral, Pierre, and Lawrence M. Murray. \"Sequential Monte Carlo with highly\n", " informative observations.\" SIAM/ASA Journal on Uncertainty Quantification 3,\n", " no. 1 (2015): 969-997. http://arxiv.org/pdf/1405.4081v2.pdf\n", "- C. E. Rasmussen and C. K. I. Williams. Gaussian Processes for Machine Learning\n", " (Adaptive Computation and Machine Learning). The MIT Press, 2005.\n", " ISBN 978-0-262-18253-9. URL http://www.gaussianprocess.org/gpml/." ] }, { "cell_type": "code", "execution_count": 3, "metadata": { "collapsed": false }, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "True posterior of last state: Distributions.Normal(μ=0.04290846886466264, σ=0.009025242808127828)\n" ] }, { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAhcAAAF6CAYAAACqW3pRAAAABmJLR0QA/wD/AP+gvaeTAAAgAElEQVR4nOzdeXwU5f0H8M/3md2cRC4BQSAb8EbYXcCzteKJCWLl2CAeiLSicthaFVq1iq1V22ptK9jWHkLVIllEfyIgHpVqaz1CdiZpKsqxuwEFpFxCzt15vr8/NtGIZMKRzUT8vl8vXi92kn2ezzzPzJdhZ3YGEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhCvI7QAipbhw/CzFWP3MS4tf+8LyouJxID2UYCxctGzRf4oLx89iovLw8vBLh9rX1GFTvbt67fh5kvWvalG7LU/l3m2D//zssmc3HP6afHWNuXhMT6/H+KGXM+c8veLpT93OI9pW8+1+yYolm9zO83VUWFiYmUe5D7BBj4aXhqNt0eZXfV4LCwszj8QarNwOIBoRTWeibzVfNKFo/HcBDgMqc9GyRf9p+j0FnH84XW3us9kL0C1ej7en1+vNANN0L+jYw2nzSOBRnq4A3VKra3MPp50JRaGnQ4Whm9oq15Gqvcep+XbfXn26ob3H9WD683q9GQDdggTarN64Pa+HO95Hag32uB1A7N+EovHfZdDjxLh/0fKSu9LVzwsvvLAHQJd0tf91pIETFLDO7RwdnYxTerT3uH7d5/Fw1/9IrcFycNEBNR1YAHz3ohWL7zvc9iaPmJxVk1NzJ8DfAGg9JenPTPqzn9XmVD8OUg+wtq9lou3h5eFfNr03dEnom6R4Zq1Rf93SpUtriguLrwfxWADZYHptD/b8YsWKFfWhUMigGjzDjJ8pxngodF60LDyzsLAw8yh0uoMJ5wBcQWT8EVrPTGh7zpKVSzY39ePQbjbV4K8EdTezfT2DziZgEwP3hpeHKz5bv+y9c0A0EsAWIl76qa7+84oVK+qd2m5pvDyG5/jiouLfMPhEBVgJtu/Y9+PWltosLix+BOCBTAgVFxYnoXQWWH1csrxkXuP7pkBxYU517jXzV82vGz16dE62zlpArO5ctHzRhweS9XDGan9aXJdRxdeCdZFheKcvXLrwfwAwoSj0ABgNi1aE7wmNDHUjxd8F0UUAmAn/tLX9l+Zj1dLc5CHvwebjVLKi5Kf75nJq/4C3ixa2+/1papMJv1Qa4xk4G6B/cQ3/VOXydcw0HsAeZvwyvCL8j8b3ZFA1/RjgkQC8IHqtzqi794UXXthzKO0d7vzuu/3tO66HtI4O87C//lrbHwFAKe4+oSj061T/qGLG3eEV4f+2NgYHO69jR47t7TWMnyutH7KVcQPAZwBUbniMHy18YeHW1uYQAMaNHneSYRuTwDwShE8IPG/R8sUvOo33wdTJ7Orc25tq8KJli95vZZvab51tcaN2kZwW6WA+P7CgB0uWH/6BBQBU51T/H8AzCXgH0PVM+v+afrbjqB1eBq5h8DGAqibgjlAolNH0c6UwlYm6L126tGZCUehhED8I4jIG5jPx+E7U6RUAtG3bNgJjPAFPM2E8mC0AyKNOz4MwEaClBMpl1q8wYarX4/3sSN2p3fr6eg8Y45n1SrAiAj0J4AwCljW9vya7+nkiGsvAbwB+g5nm5Km8n7XWdssjxosBvR3AAgYu9JDxr8kjJmcdSF42+A0Au8FYC4V3oWEAfONnTRNfD8b4uk57zwSALDvrG2AUfsqfxg8k6+GO1b6c2ktkJF4E6FvaTv4GACaMGn8FA7cDagUAkIG/sqJJBF5EzAuIMcZDxrPN229pbr40Tvvh1P6BrKvTdr8/TW2SRhiMzkx4GeAbKQcVzHQNES8Gw0MKLxYWFmYCAFVjCcDjmPEwEf8czBdkJTNfnTNnjjqU9g53flsb10NaR4d52F9/TvvjZzlBv2fiPIAWADiV6PN1aG0/OJh5Ncg4ioFrbKWWgfERiB4FeIhtJ//ZVOec5hAAGbZ6FkBngOYA/CGDnr2icFywpfU/2Dr5xRrsnKelOtsRyScXHQgTxgA0CMBuMI8IhUJGOBy2D6fNUGHoXAIuZsaIksb/iRQXjo+C6KF9f1exelKTfa/ayxcAWDF5xOSsGlRfrjTfFBoZOo6B7xHx1YuWLX4GAMYWjn3VQ8bG4lHFIz6p/uTN1Ergg5IV4XEA+IpLx4/QGpckDX3ys0ufXZPqu/gOEH9WaFprtw51panfpAUlK0p+DAChS0IWKbwZKgodYyg+SWuM1ITh4WXh1QBQXFT8KTNPbq3tkmUlr+930IhuL1kWXgAAodGh58jGhuqc6isAzG+1zRdLngsVhe5QgLloWcnK4qLiWoBnjRkzpnun6k519agbCsa/tVYjAKwi5guY6NVOulM/NpyzHu5YhZeHt3xh22htXZ4reX1C0fjrGbS0uLD4NWb+OTF+tmjForcBgIAdDNyzaPni1Y3bWpIITzW13zj/+52b8LLwbc3Hab/T0Er7TusKxokHut1/qV/G/EUrwvcAwISi0HEMXMXAWSXLFm8JXRJ6mRhrOiPnlFBh6CgAo5gwPLw8tX5jC8f+00PGxsrSym9D4dWDam9kaM/hzm/4xfBzrY3rwWQCEHGah337c5rzffr/86Lli+8CgCtGjS/XTG9MvGxir2R9Ms9pDFizPpR5ZWBeeEXJgwAwduTYlz2GEUM1JoQKQ1VOc3hF0biPNXAKg68KrwibAJaGCkNrNBl5+1v/Q6mTl112WV5TzlBh6FzHbQpYuu/7ndbbTXJw0bEMAtP1TLyRCCvVXtwN4J4DfXPxJcVFTHz7ZwsYNxLxOQDt6fpJ17eaFpMyVjLrL+2Mz7z0TKy4KPQGK4wHsKI2d28RmFh3oiVUQ+PArFirAaGi0A+avW0nMwcBvJlqmxeicYPXrIYD/FHTgQUAsOKVxPjs4II8dGYr7ZYCALF+s+kHhjbiWtlQhupj2/o0IuwNLwuXNf288RTEvOJRxVe30vZ+Dy6UoV5p+nt4aThaXBSKEvNJB5j3C21+UvPJWz1zenzqrfOeU+epqybgQ1K8kBnjAICJLgDhjwfS7uGOFYAvHFwcSJ+Lli9+cUJR6M9M/GcA722t3fbZx+yLlocnTbh0gr94VPFtYO0DMBrNPg21NbU4N/sb93211r7TurJtH/B2/2X89md/A/8PoMhnB2aMbQCgldEdqTHaPei0QZHwsjAAYMmKJZuKi0JVYAwGUgcXB9oeKQxqy/ltk3XEgc1DkwOdc2Z66/M3eTZC2dCsjyUPneI0BkSccyjzqlitavr7kpVLNhcXhT5QwCmacAwc5rAuw349o0FFifFq8ajxS8DqVW7gv5W8WrJ7f/0cwD71pTr5xQZwulMeNB5ctPj+DkROi3QgBDxasqLkL+Hl4VcAfpQJdxYXFn/jQN/PHt6pFP7T9Ade1ALoRsBHj69+PNH0e1rrrS1n4AVgunzEiBEeZnUFiBeFw+FaMHoAsElxb0UoaPpDhKcV8/uftc206/PWdFcAO5q3byftj78Y+oDbrWn6e4PR8NlORaCjkSqqX97RDrDtL48BfeHTIgb2gBQdSpurVq1KAvQySH+LiEeA6XWt6TUAZ44dObY3AUFb28sPqN3DHKtDHR8m3to4MNtS65NSXBR6kbVeCa1PZVAcRI/tM44tz80BaK39Vtb1oLb7fdrcs8+i6v39XuP6Vc2ZM0d/cTlYaRgH216bz6+DA86EA5uHJgc65wze/1e9Wx+DQ5pXYvrSgRcT69bm8Pnnn9/FNRhCxLczU0+A/0QZ+GjCpRP8h5g/lfkLdbJ5vwe8Te33/R2JfHLRkTA+OxrOqek0uyan+iIifjp0YcgffjW83yPl5sIvhv8N4N/Nl00YNT7KjBNDhaEe4RXhbQBAROe1tO/XehoWZyUz5/bo1GM0NI9SwEUAwMTrieFRSf37hSufrWz8dQqNCk1RMD7c/+pQnIBBV466suvflv1tJwB4lGd4874Ppd19RMEYOHH0xKM/v+hw/LcBGqeJFx9K27ZtnwJgKwBcfvnlXagBp4Lxm0PNS9ArmGkGgFo2+OHwi+H/FheFdniU50cA/3fJiiWbQqNCrbbbBmP1BQfS3hVF487QTLNB/FMw3VFcWHx9yYqSP15RNO4MDRQp1sOeWfFsBGi6iO0LXbQ4N4uWhyc5ZTvA9ltEhIPa7g9RFMCg0MhQt/DK8A4gdQEhA/lM+OBgG2vr+W0LhzAPhzznQOtjQGQfdyjzaiv7eAAxIHU/GwCnksZDTMiBwxxOvGTckKQHAxe9uPgJAE+kLiatfgFazwJw1cHmb2390cbblJvk4KKDmr9qft2E0ROuZlu/TV76PYCJzX7c/4pLx49o/vtsG3VN58Kb00xhItxPoFsA3DmhaMLxDH1zS/2+8MILeyYUhZ6Hxm8BbHxm2eK3AKBO1b2abWdttJVxV+jS0G1e9tYkOXk7wDMTOnHC/toyDONFnbQfSqLh56Gi0N0AehD44eZl4FDa3aePxbad/Llt278Yc/GYH6oM1Zk1/RygN+tU7SG1zcwPjLl4zGWezp7tVMM/ALCrzlsXPtC8BNSAcGzqr2ANWkGExwHYSKDxWwH8GghTCXjkQNs93LHaV2vtjR49Okfb6q8Ani9Ztvju4qJiDxT/KjQ69CpDJWBrYngyASBUGDoFxLcBoDlz5qg5c+Zop7nZ3zh9YQ4Mj2P7ZWVlcHKw2/2hMDzGc7ad/KVS+B5Spy/J4/HcCuZP6oy65wlktNZGc201v07jerBam4c5c+bo5v21NuetaW0MDDIOaV4JuCd0YejdZF7S420w7gNjpyZ62eMx2GkOO9kZBaTxdKgwVBheEX6juke1pmp0Zv58fZqv/2HXs1a2qQMZw45CTot0YIuWLioD0xwQXzGhcPxnR/0MTNCaXm/+h0kv3F8b4eXhLcR8NcAziotCuxj6P2B+DkBtix0zLwDQl4H5TYuWLl1awwrjQXwmaUSTnPiIwBMAGt/8K6XNPfPCMx8zMA5M5xPwMQEvg/gxACCbag+13eYWLl34PwbGAzza6/FUGVq9R8DqhoyG2w+1bQL93evxVFE1PoWmaUrxuKavpR1Im8x4hoGriotCLzfNARjlAN5v+t8Ig14DkEnEyw+03cMdq3211l6Ozvo5gJ5J254JAHt4z71gbISNBScPO9kE8wIm/Y/iotBWIrwIptsAfPLfdyr/2drc7G+cmktt+87tOzmk7f4gLXxh4VZinsiEmcVFoe3FRaGd0DyRFS5funRpTestfFFbza/TuB6sA5mH5v21NuetaW0MDnVeCVhDGdjirfdsJKYRTJgUXh7e0tocLlz5bCUBi4iwqrgwtImqsQtAJ8Nj/K6p7ebrf9j1rI23KTfJ7b+/JkIjQp2MHD1IUeaHTacoDsWIESM8vTr1GoQkMjpv62w2P/f5pT5DIQPbkB1eFd57zcXX5D758pPVoaJQgIA3SpaHO6PZ/6oOpt2W+tJ1+niVpTaEw+GGQ838WXsjQp2Qg4JtNdveb36dwYG2GRoR6gQA4VXhvQezHgeS9XDHqi3bG3PxmJ4e8hwVXhleB6TuQfBp5qddmxdSp7lpbZwOpH0nbbXdO5k8YnJWXc6eU5k81TpHf3i43/Bqi/k91O2vJa3Nw779Oc35gTiQ/etA5jV0SehEUlhjkx6YMBLbcuu9/ZudrvhMa3MYGh0qIE0+w7a3J/NU5Zd+vs/6H+4ctvU25QY5uBBpEyoM9SDCeiKevGjZ4uemDpvq2dVrx1MAdS1ZHr7Y7XxCiCNb84OLI+m5HV8FclpEpE14RXgbM2Yz0x+Ki0Kbd/XauRuggYlk8mq3swkhhBDiq40mjhw3aOLoiUe7HUQI8fUxZ84cdVXhVUc13m1TCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIY5IFRVVAysq1vVzO4dwFolER7idQTizrC25prnhdLdzCGerV8dOrqyMHuN2DuGsPWue3FikBaFQyJg+afoJoVDooJ5sCADJpB6ptfGNdOQSbUnd63YC4Uzr2h7MNN3tHMKZUhibSMDvdg7RmvareXJwsR+3fOeWbr1ye/6HFO7tldvz3WnTpnU6mPcrxWW2ze+nK59oG0T6b25nEM4MQ+0mUsvcziGcEeEtQEfdziGcSc1z2YxrZ9w5c/K06wFg5nXT75157U2j3M4khBBCfFV43A7Q3m4J3ZKdyE1cOnf+3HDTsllTpuRV29mFRNTgrfauTFDDKZroVQDQGpVENBTAAf/vybKiZ9o21w4dOsBKwyqINmKa0RsCgYI/uJ1DtKy8PN5Va1wUCOSXuJ1FtMw0Y+cB9keBwMAP3c4iWtaeNe9rdVpk2nXT+iVy6x8A8S1Ny6ZOneqt0dlvAggBuCiRU/8CAw2KVB0AgJEAqdqD6UdrGmoYdHKbhhdtjlld6XYG4cy2dWdmLZ8cdnDMOBtQBW7nEM7as+Z9rT65UJoeA9HRAHPTsoxExgQQx+fNnxcCgBmTp0dA+IfWeggAC4TBpPXbB9OPx6NWAg0NbRxftDl9j9sJhDOlsrcx18xzO4dwpjWWZGVhp9s5RGvar+Z9rQ4u5i6YN3r6pOlnkYGHP1vIGATQ56cvmN43iLdowk3Tr51+GSnO7Z7f8/592yot/fjojIzEUQ0N3i3Dh/epsaxNfYnsjDVr+sdPPRUbKiqqCioq1vUbPPi4jW+9tTG7Uyfdmzmx1+8/7pNIJNrFMFS3mpqG/5155vGfVlZGj7FtlZOd3fDR8ccfX19WtiHfMJTy+33R0tJSb0ZGj37JJNUNHdr/48rKTzrZdm1Pw1C7Bg3qt6MpB3PmVr//mGrT3HisUjpzzZr+8VAIuqKiqqChIZkcPnxgVTQazdqzR/XZN0dDQ3L78OEDd1dUbOjFbOQ2z+HxGMaQIfkbXn+dPd27V/XfN0dtrbH7jDP6bn/nnU3ds7Ptzvvm2L69f9V551GyvDw+YN8cRHb14MEDtpaWru+ckeHpvm+OvDz9cUFBQV1p6fr+GRkeT/McWqv6QKDfR5a1JZeovte+OQwj+5NBg3ruLSur6uPxcFbzHMmkbQ8dOiDeqVPy3+Xl8QH75rBtvSMYLNhlWet6Enk77Ztj8OD+0XAY6qSTqvKbcpSWfpyTkZE4pqHB++nw4X3+V1m5sZtt6y775mho2LZx+PDhCcuKFdi21kOHDoivXbs2s7Y241jD0DWDBhVsefvttUfl5GQcvW+OvXvV5rPP7ldbUbGuH7PX2zwHs9Hg9/fd1FKOjIzabSeddNKe0tJY74wMym6eQ2vmYLAgVllZmWHbnfrumwPAziFD8neWla3t4fFk5DXlaNrmBw/uHwOAiooqX1OOpm0+mWzYM3To8dvKy+NdAXRtytG0zRvG3k2DBg1qiESiPgAIBgtiTdt8QwPX+v3HbF6zZs375eXxAU05Wtr3GnNwRUVVAVEi0Xzfa8rR0r7XlGN/+15DA9cOH+7b3NK+15Rjf/teU46W9r2mHPvb95pytLTvNeVoad9bs6Z/vLiY7P3te005Wtr3mnLsb99rytF83xs2rO/777yzqXt5eXxAU4797XtNOZq2+ZZqQEv73pAh+RtKSthovu/tWwNa2veacrT1vnegNWB/+15Tjpb2vaYc++57zXMAoOb73r45vrjvFaxqyrG/fa8t/739Wp0W2S/mXgSONb0k4s0Mzjo6v8c3vXbi5rlPPFY0Z86c5L5vM4yGkVrr2RkZ9sBUM8nJWuvZJ564Nes//9lYzKwftW3PVADIzU3011rPZjZGpd5LZ6Xe6wkCQEMDjdVaz969W/UCAKXUTAC3pf7euYvWerZSeiIAJJPVg7TWsxMJ+1wA8HrrL9Raz9a69vhUfvtarfXs/v3X5YbDUFrr2YZh3AAAO3Z4+mqtZwOe0anfxRla69lerzEMAGybxmitZ+/d6+mdyqmma43bASAn56POWuvZRHw1ANTX152ktZ6dlWWfBwCZmfYFqbbrTkyNUPIarfXsLl1inZiZtNazPR41DQD27FF9tNazbdv4dmodPael1tFzWiqH8e1UDuNYAPB41LTU+DF16RLrlOoneU2qn7oTtdazMzPtCwAgK8s+T2s9u6Gh7uTUOuqrtNazc3I+6gwAWuN2w1DTUzk8C1JjR2NTY2kMS60jzki17RmttZ69a5fulxoP4wat9exwGKp//3W5qUz25NQ61B2fGsv6CwEgkbDP1VrPTiarB6V+riem1rFzl8ZN6LbGeUZdHfVMZaZxAJCR4Qmm5o3OSm1bxiit9ezc3ET/1Ph4pmqtZ69evdpz4olbs1I5ktel3msPTL23YWRq20qeo7WeXVOTdWpqHXhCary6dUu1jVuUUt9LZT6qR2rbUiEAyM7O9Kfmic9OvTejMDXnuiA1lonvaq1nx2KxjH//e1Nmal7s76TmQRek1jfjksZ+vqm1nl1Xl+1PbceqWGs9u74+7+jG8fg+ERpPVx59dGodeEJZWVWfurqs36TWkc5JbQ+JS1JjnRyQatuekhqPzVnr1q3LSGX2XA8AnTol81Nz6mk8tUJnp7YHI5BaZxqf2gYym/a97zHj1tT20KVrartF0753amqekt9KjUfDxal1rDu+McdkrfXsIUM+yFm1CkaqXc8NALB7t26sAZ5LU/3gzNRYeoelxjJVA/bsMY5JbbdqRtO+l5e3uUtqHfRVqcx7T2msASNS+17iwlTbdSek3pucpLWePWDAhk7M3FgDPDcCwN69xrGpXN5vp9bBc3pqLD3DU9sWXZ7aXrx9Gve96VpjFgB0716Vl1pfvibVT/1JjTXwfNOMTcnISN6Yyll/UurnfLXWenb37lV5jfveLI8nte/V1HgbawBdnhprz/DGnI33NPFcprWevWOHp29jzhsba4AaMGBDp1Q/yUmpdmuPT23TiaZ9b0RqfPaeknqvvlJrPTsvb3OXxhy3E6kZqXlRvRr3vbGp7dY7LLWOODM1p55LU9tLQ7/Uez03aK1nr1oFY8iQD3Iat/nrUvNQd1zjvndx4773rcYacGpqLDExNdZdujbuE7c27XvJZOeeqcw0HgBycoxAqm06OzWWnlGNdc6XmifP9Vrr2evWrctYvXpzYw2wp6TGMjkgtd0mGvc9OkdrPbu2NntIJBL7o9dLxanXud1TbeMWAN+HODzTJ00/a8Z109/67PXk6d+fMXnaj5u9fuLmSdPOO5w+IpHYNMuKXnE4bYj0i0Ti/3A7g3AWiUR9kUh0gds5hLNIJHanaUZHup1DOGvPmve1/+SCmN4D6HwAuOnKm7oCGGGDyg+nzfr6xFMeT+6LbRJQpA0zTXQ7g3C2a5dvk1L27W7nEK3heTU1njfcTiGctWfN+1pdc7E/W2u2vt0rt8fWmZOnv8tAf2Lc/+hf520/nDbPPPP4T9sqn0ifoUP7f+x2BuHsvPMoCeATt3MIZ8FgwS63M4jWSc1zwYwpM/rMmjIlry3aktMiXw1yWqTjk9MiXw1yWuSroT1r3tf+k4smc/8yt82O6IjoU2ba21btifQgYjmK7+CU8iaZk9vcziFatZPIqHE7hHAmNU8IIYQQQnyurGxtj8rKjd3cziGcWdaGE93OIJxVVlZmNH0PX3RcpaWx3m+/vfYot3MIZ+1Z87723xZJByJvKJlMXux2DuFMa+NxtzMIZw0NOX0AtNtjosWhMQxMycrynOV2DuGsPWueXHORBkR6ndZKrrno4Ij0KrczCGfMyWqlvO+5nUM4I8L7WvMWt3MIZ1LzhBBCCCHE50wzFly9OiZPRe3gIpGYPBW1gystXd85EonJU1E7uLKy2FmWFZOnonZw7Vnz5JqLNGDGWR4P+93OIVpDN7idQDgzDNUV4GK3cwhnRDifmU9wO4doTfvVPLnmIg2I8O9kkurcziFaw39wO4FwZtt6p2EYJW7nEM6Y8XelSK656PCk5gkhhBBCiCamueHiSCR2tts5hDPTjMpXHDu4srK1PUwzNsPtHMKZacbGlpVtkFPBHVx71jy55iINmNVxSnF/t3MIZ8xqhNsZhDMiTy4zn+Z2DuGMGScrRce4nUM4a8+aJ9dcpIFSnheUaki4nUM4Y6Zb3c4gnNXWerbm5OhfuJ1DOFMKf6ur88rToDs4qXlCCCGEEOJzphm/2jTjl7qdQzgzzdhCtzMIZ5a1qa9pxh5yO4dwFonEpllW1Tlu5xDO2rPmyWmRNGDmo5TipNs5hDNm6uN2BuFM64QHQA+3c4hWdWW2c9wOIZxJzRNCCCGEEJ9bu3ZtZmVlZYbbOYSzNWvW5LmdQThjZmVZW3LdziGcRaPRrNLSUq/bOYSz9qx58lXUNNi71/udZDJnrNs5hLPa2uwX3c4gnJlmrL/WtY+5nUM427WLbvV4up/vdg7hrD1rnhxcpIFS9Akz73A7h2gNf+h2AuGM2WggQsztHMIZETZrTfJV1A5Pap4QQgghhGhSUVE1sKJiXT+3cwhnkUh0hNsZhDPL2pJrmhtOdzuHcLZ6dezkysqo3KGzg2vPmienRdIgmdQjtTa+4XYO0Rolzxbp4LSu7cFM093OIZwphbGJBOTZIh1e+9U8uc9FGjBzpdZGjds5hDOleKnbGYSzhgbvnszMxCq3c4hWRbRWG90OIZxJzRNCCCGEEJ+zrOiZ8vjhjs80oze4nUE4Ky+PdzXNeLHbOYQz04ydZ5rrT3A7h3DWnjVPrrlIA61pqGHQyW7nEM6Y1ZVuZxDObFt3Ztaj3M4hnDHjbEAVuJ1DOGvPmifXXKSBx6NWAg0NbucQrdH3uJ1AOFMqextzzTy3cwhnWmNJVhZ2up1DtEZqnhBCCCGEaGJZ8dFyD4WOzzRjv3I7g3BWUbGhVyQSn+V2DuEsEoldWVZWNdztHMJZe1LpLAoAACAASURBVNY8OS2SBlpzP6Ugt//u4JhpmNsZhLNkkrIBPcjtHKJVBUrZ290OIZy1Z82Tg4s0qK9PPJWX11m7nUM4Y6aJbmcQznbt8m3q1m397W7nEK3heTU1nnq3UwhnUvOEEEIIIcTnIpHYNMuKXuF2DuEsEon/w+0MwlkkEvVFItEFbucQziKR2J2mGR3pdg7hrD1rntznIg2UQj0zyVdROzgi3uN2BuHMtrUmUnIr/Q6OCLVElHA7h3AmNU8IIYQQQnyurGxtj8rKjd3cziGcWdaGE93OIJxVVlZmRCJRn9s5hLPS0ljvt99ee5TbOYSz9qx5clokDYi8oWQyebHbOYQzrY3H3c4gnDU05PQB0G6PiRaHxjAwJSvLc5bbOYSz9qx58lXUNCDS67RWe93OIZwR6VVuZxDOmJPVSnnfczuHcEaE97XmLW7nEM6k5gkhhBBCiM+VlcVPqaioGuh2DuHMsuKj3c4gnK1ZsybPNGPnuZ1DODPNWNCyNvV1O4dw1p41T665SAMiHqG1fZrbOYQzrXGb2xmEs9razO7MPNntHMIZM4qYE3Kb9g6uPWueXHORBkT4dzJJdW7nEK3hP7idQDizbb3TMIwSt3MIZ8z4u1Ik11x0eFLzhBBCCCFEE9PccHEkEjvb7RzCmWlG5SuOHVxZ2doephmb4XYO4cw0Y2PLyjb43c4hnLVnzZNrLtKAWR2nFPd3O4dwxqxGuJ1BOCPy5DKzXL/UwTHjZKXoGLdzCGftWfPkmos0UMrzglINcp/9Do6ZbnU7g3BWW+vZmpOjf+F2DuFMKfytrs77qds5hDOpeR3E96+e2nvmzJmZbucQQgghvkrktEgLvjflphOTHu+rtIcO+vSGacavNs34penIJdqOacYWup1BOLOsTX1NM/aQ2zmEs0gkNs2yqs5xO4dw1p41T06L7MeN19zYM6nVLQByDuX9zHyUUpxs41iijTFTH7czCGdaJzwAeridQ7SqK7N9SPVStB+peR3EjMnTnrt58s3Hu51DCCGE+Cr5WnxycUvoluxEbuLSufPnhpuWzZoyJa/azi4kogZvtXdlslPDD5lxKQNPzps/79eH09/atWszGxoaeNCgQQ2Hn16ky5o1a/JOOumkPW7nEC1jZlVevjXb7z+m2u0somXRaDRr+/bt9vDhw+VC9g6sPWveEX9wMe26af0SXH8riE4HEAaAqVOnemsavG8CWAvgk0RO/fS5Tzx2EYB72qLPvXu931HKswPAM23RnkiP2trsFwGc63YO0TLTjPVH6pHr17qdRbRs1y661ePpXgpgpdtZRMvas+Yd8Rd0Kk2PAeqM5ssyEhkTQIjPWzAvNHf+3OkgOnrGNTMK2qxPRZ8w8462ak+kC3/odgLhjNloIELM7RzCGRE2a03yVdQOr/1q3hH/ycXcBfNGT580/Swy8PBnCxmDALI+f03vk0cPAxD9wnvnPzampXYtK/ojZhpHpG7x+/u/GYnE/kaEE5JJ+4KcnIal1dXef5lm7KZAwDemrKxquFL69wAvCwQK7jHN6HUATSeiB/3+/MWmGfsVgG8pxVcPGVKwxrJirwDI9ft9Z1vWpr7MyeeJUOb3+6aWl8cKtcZPAcwPBHxzI5H4LCIuZubbgsGCVaYZexLAyV6vcfEnn/T9tGvX+NsAfxwIFFxmmrEggD8S4SW/33eXZcUnMfPNRPyQ31/wTCQS+yURzgPo2kAgv9I04y8BulsgUHB6aWmst8eDpQCsQMD3nUgkfhERP8DMTwWDBb82zditACYS0Wy/P/8104w/AfBgomTRkCEDt1lW/D0ibPX7faPKyzcO1tp+AuBXA4GCH0YisauIcAszHgkGfU+bZvwBgC8yDD1l8OAB5ZYVW8aMXn5//mnl5et7MHuWE9F//P78yZYVO58ZvwCwMBDwPRyJRL9PRFcD+o5AYMDLphn7E4BAMonRw4f7Nptm9F1A7QgE8i9hpkdMM1bKjNeDQd/tlhW9gpluI6Lf+v35f7Ws2H3MuATA9YGAL2Ka0RcA6rNzZ/6ZPXtuOiqRsF8G8H4g4LumrCx6rlL0MDOVBIP5v2i8q+RkpfDjIUN8K0wz+geAhhF5Lvf7+26yrNhbAKr9ft9FlrXhRGb1NBG96ffn32JZ8fHM/EOA5wUCBU+k7qpHo5jVTcFg//dMM/YcgH65uYlv7N6tsjwe4zUAawMB30TTjH8T4F8T8bN+f8EDlhW9iZm+Q0T3+v35SyOR2GNEOF1rPW7o0AFx04y9CSARCPjOX71643GGYT8D4K1AwHezacbGALiTCL/3+31/Ms343QBfRsQz/P6CtyORWJgIBclkxreysz0qkahZRYQNfr+vuKwsdpZSeBTA84GA7z7TjE8FeCoz3xcMFjxvmrFHAZxlGGrC4MH911tWbBUzq0Cg4FuWFStgRpgI7/j9/adbVtVq04yVAvynQKDg95FI7E4ijGGm7wWD+f9qvAL++Lq6xPlaZyVycuw3AcQDAd+4SKTqNCL9O2Z+MRgsmGNZse8w4yYiPOD3+561rPgjzHwOYF8ZCAz80DRjrwHICgR836ioWNfPtj3PAVQaCOTfGInERhHhXmY8EQz65plm9IcAjVdK/2DIkAFvmGbsKQAnKUUXGcbe6kQi9y0i2uT3519uWbGhzHgcwIpAwPfjSCQ6mYhmAPSLQCC/pPEbMSNsG9cMG+Z73zSjLwN0VCDgO7OsrKqPUvoFZkSCQd/1lhW9hJnuI6K/+v35v7Ws2G3MuALA7YGA7/VIJLqAiAYlkxmXDBvWe4dlxd8FsDkQ8I0uK9vgV0r9GaBXAoH8H5lm/GqAv0/ED/v9BQtNM/pzgC5QyrhuyJB+FZFIbLlS6OH3+06rqNjQy7bVMoAqAoH868rK4hcqxQ8C/HQg4HvENKO3mGbsUWb6UTCY/4ppxv4CYIhh6FGDBw/Yalmx97TGtmDQV2RZVacy6/lE+Lvf75tlWdGJqXsw0K8DgfynIpH4/UR8MTN/NxgsME0zthRAb78///TVqzd383gaXgLw30DAN6nxqbm/JMIzfr/vIcuK38zMk4j4Lr+/4KVIJPZHIgS1VpcNHdr/Y9OMvQ3wp4FAwcWrV8dONgw8CeAfgYDvVtOMFwM8i5nnBoMF800z9lMAhUrRDUOG5K+2rPjzzNzX660+27Y75WrNrxDhA7/fd5VlVZ3DrB8BeHEgUPBgJBKbToTrmHFPMOhbZprx3wF8mmEkxwwefNxG04z9C0BdIOC7IBJZfzyRsZCZ/xkMFnzfsmLjmPEjIvzO7/f9ORKJziGiSwE9LRAY8K5pxp4FkF9TY5yjVJ03K8v7d4DWBQL5V0QisbOJ8FtmPBcM+n5mmtEbAfoukfqJ39//etOMzQVwJhFCfr8vaprRN4hI+/2+EQf1j6sApk+aftaM66a/1fR6xrXT/zJz8vQpTa9nTp7+8MzrprXZLYYrKqoGVlSs69dW7Yn0iESiI9zOIJxZ1pZc09xwuts5hLPVq2MnV1ZG5Q6dHVx71rwj/rTI/jChnMHHNr3WQDeyUdlW7SeTeqTWxjfaqj2RLkqeLdLBaV3bg5mmu51DOFMKYxMJyLNFOrz2q3lfy4MLYnoPoPMB4KYrb+oKYIQNKm+r9pm5UmtjfVu1J9oWMytmHtyr11EfMvNwZs5yO5PYv4YG7x4irHI7h2hVRGu10e0QwplSvLS9+qL26shNTddczH1i3tkAEAqFjF65PZ4m0AAG+hPj/kcXzPut2zlF+jGzr8bGHfUag6rqkX1sJhpyFDblGHgAgElEttsZhRBCfIXNmDKjz6wpU/Laul3Lip4pjx/ueJjZW5Pk8DOf8M6LTU4OK2X7m2Wc/M1Grt2Z5H9WM/d2O6P4ovLyeNfURXaiIzPN2Hmmuf4Et3MIZ6YZvaG9+jrivy3iZO5f5n6cjna1pqGGgR0ArFZ/WbSnwF4bZ/2yCkdZe0FJTn1yF6+HOq8b/Cdm4ptovBeK6BhsW3cGMApAidtZRMuYcTaRKgUgX+/uwJjVlQD+0B59fa0PLtKFmVYpperdziG+5LiqenjX1X5+YAEAG+tA71fDGJiFE90MJ74sO7t+e11d1ny3cwhnRFhO5N3mdg7hTCnIQwCFaGvM/O11Nbyl65us8Tpz8z8rt3N1jc23uJ1RCCGOBF/Lb4ukm2XFR8s9FDqk93pmoCbUE5ylwE0LC7sBp+SCPQr/dDOc+LKKig29IpH4LLdzCGeRSOzKsrKq4W7nEM4ab9jYLuS0SBpozf2Ugtz+u+PZkmNgzs3H4r4Lu+DoDXXI7J2B+kG5SBztxYO7/4c1bgcUX5RMUjagB7mdQ7SqQCl7u9shhDNmGtZefcnBRRrU1yeeysvrrN3OIb6IiDQzPzcwE7G+mThrV4M+Ic+jqrIMvJtUKD36aOx1O6P4ol27fJu6dVt/u9s5RGt4Xk2NR64z6+CYaaLbGYQQQgghDkm6b6JFxaOKRzBzUDG/DxgfJnSiZsnKJZvT3K+rGh+OtDMQ8C1xO4tomWnGXgwEfJe6nUO0rLR0fX+Px/OjQCD/JreziJaZZuwHRDD9ft/f3c4iWtaeNS9tF3QWFhZmhopC/wDz3wm4C4qKtNIjPB7jw+JRxSPT1a8QQggh3JW2Ty6KC4t/zMTjbbZHGWQUK0LBomXhm4uLin9G4Amdt3Y96fHVjyfS1b8QQggh3JG+r6ISXwDgkSUrlmxqtpQ95PklA7139t4ZTFvfLisrW9ujsnJjN7dzCGeWtUFumtXBVVZWZkQiUZ/bOYSz0tJY77ffXnuU2zmEs/aseem8z0UtER+/78IEJXIAeBSO3DtYEnlDyWTyYrdzCGdaG4+7nUE4a2jI6QOg3R4TLQ6NYWBKVpbnLLdzCGftWfPS91VUpiUAHpwwanwFa8oCMSaOHDfItukuEH+0de/WyrT17TKlaKPWvMftHMIZEa92O4Nw5vFwbTJpHLG14ggS1dqQ+1x0cO1Z89L6bZEJRaGHGfgeUp+Q2EgdzFSxQij8YvjddPYthBBCCHek9fbfi5aHb2UbJ4FoEkCzmPDtWqPu5CP9wKKsLH5KRUXVQLdzCGeWFR/tdgbhbM2aNXmmGTvP7RzCmWnGgpa1qa/bOYSz9qx5af0qavGo4odgYGHJspKnSpaXPKIY47Pt7N9cfvnlXdLVb0dAxCO0tk9zO4dwpjVuczuDcFZbm9mdmSe7nUM4Y0YRc0Ju097BtWfNS9vBxVEq9ztgnkSg+Z8v5WcBPiezwTsvXf12BET4dzJJlts5RGv4D24nEM5sW+8EqMTtHMIZM/5ORB+6nUO0pv1qXvruc1EUWgHwGyXLFz/QfPm4S8cdb2j1HueiZzgcbkhX/0IIIYRwRzqvuTiGQV+6zbdu0LsBZOk6PSCNfbuqrCx6biRSJadFOjjLislpkQ7unXc2dTfN6HVu5xDOIpFYUVlZ/BS3cwhn7Vnz0ndwwbxSAXeERoWGofETkomjJx7t9XgfBLBHZam1aevbZUQ0SClbLujs4LQmuaCzg8vISOQxY4TbOUSrgkrpfm6HEM7as+al7T4XObWd5tRkVwcIKC0eFdoLYIttJwcC2MPA+HA4bKerb7cp5XlBqQa5tXkHx0y3up1BOKut9WzNydG/cDuHcKYU/lZX5/3U7RzCWXvWvHQ/FRVXXDp+hK3pNAIdDSBqGMbihUsX/i/d/QohhBDCHWm9zwUA2JrWKEO9RgYtIoPe1dD9J4yeMDTd/brJNONXm2ZcHuXdwZlmbKHbGYQzy9rU1zRjD7mdQziLRGLTLKvqHLdzCGftWfPSdlokVBQaTMAyAP3Y1vv7lbR/auIWZj5KKU66nUM4Y6Y+bmcQzrROeAD0cDuHaFVXZjvH7RDCWXvWvLQdXBBwDwBijXMMbcQbjAZOV18dTadOiT83NHx91verKju7Vj5d6uACAV9VefnWaW7nEM66dOGHt2/ffsReR3ekaM+al7aDCwb6EujB8Esl/0xXH0IIIYToeNJ2zQWB3wLjuHS135Ht3ev9TjKZM9btHMJZbW32i25nEM5MM9Zf69rH3M4hnO3aRbd6PN3PdzuHcNaeNS99p0XI+CNDrywuKt5LrN/UTDXNfx5+KXzEfqKhFH3CrOVrWR0ey+2KOzhmo0EpO+Z2DuGMCJu1Jql5HV771bz03f57VCgMxviWfl6yPHzEXtAphBBCfJ2l7bQI52BSnaf+qJb+pKvfjsCyNvWtqNjQy+0cwllZWdVwtzMIZ2+9tTHbNOPytM0OzrJiBe+8s6m72zmEs/aseWk7LRIOh2vT1XZHp3XyMoB2AHjG7SyiZUT8MIBz3c4hWpadnezFjFkArnU7i2iZ1rgyMzNRCmCl21lEy9qz5qXt4KIRFY8qHsHMQcX8PmB8mNCJmiUrl3zpgWZHEmau1Nqoaf03hZuU4qVuZxDOGhq8ezIzE6vcziFaFdFabXQ7hHDWnjUvbdc9FBYWZnaiTq8QcA6AnUR4WgMmAb8GaHzJshI5whVCCCGOQGm75iIPebMAdE6y3Y+B+wAgvCz8FzA9SsyPTR021Zuuvt1mWdEzy8o2+N3OIZyZZvQGtzMIZ+Xl8a6mGS92O4dwZpqx80xz/Qlu5xDO2rPmpe/ZIsQXAHhkyYolm5otZQ95fslA7529dwbT1rfLtKahhkEnu51DOGNWV7qdQTizbd2ZWY9yO4dwxoyzAVXgdg7hrD1rXjqvuagl4uP3XZigRA4xPAqqPo19u4qZVil15K7fkUIpyAOxOrjs7PrtdXVZ893OIZwRYTmRd5vbOYSz9qx56bvPRWHx9SB+kIins1YDSHFvldS/t5VxF4jP/KRm2/GrVq2Sh3sJIYQQR5i0nRYpWVHyRwLmM9NTIL6PGTfahvoPiM9mhQlH8oGFZcVHRyLREW7nEM5MM/YrtzMIZxUVG3pFIvFZbucQziKR2JVy35iOrz1rXvquuQCwaHn4VrZxEogmATSLCd+uNepODr8Yfjed/bpNa+6nFI5xO4dwxkzD3M4gnCWTlA1ouYlWx1eglC030erg2rPmpe2ai7Ejx/b2erxdGl+uBgEEQo7Oyb/ikitqn3npmTiAI/Kx5PX1iafy8jprt3MIZ8w00e0MwtmuXb5N3bqtv93tHKI1PK+mxiPXmXVw7Vnz0nbNxYTC0B+YMNXhV/YAeDuRTF793MvPfZKuHEIIIYRoX+k7LaLUowC2gOlOJgxP2nYfMF1GwFoivoWUOgfgnV6PZ3naMrjEsmLfMc2YPHK9gzPNmDxyvYMrLV3f3zTjv3M7h3BmmrEfWFZMHrnewbVnzUvbaRFm+3oCLVu0ouT+ZouXhopCMWJ6uWTZol9PHTb16l29du6eOHri0QuXLvxfurIIIYQQov2k76uoRaG/A7S0ZHnJI82Xh4pCxxDwcWYyK+/Jl5+sKS4KbTZsfcHClc9WpiuLEEIIIdpP2k6LEONNAt80bvS4k5qWXTnqyq5E/BMAG558+cnqCaPGjwFgLFz57H/TlcMNZWVre1RWbuzmdg7hzLI2nOh2BuGssrIyIxKJ+tzOIZyVlsZ6v/322qPcziGctWfNS9vBxafYez8D6w1bvV9cFNpUXBQqT3Lif2C6gpinXXPxNbnMNJ8Zd+MI+9YIkTeUTCYvdjuHcKa18bjbGYSzhoacPgDudTuHcGYYmJKV5TnL7RzCWXvWvLQdXKxYsaK+ZHm4UCk+D0S/BvFSJtykPMZJi1Ysfrmuc10dMwaGV4Q75MVaU6dOzZl+7fRBOIRTR0rRRq2xJQ2xRBsi4tVuZxDOPB6uBZScMu34olob290OIZy1Z81L31dRR4UeBbO1aPniP6Wrj3SZee3MgUz6ZYBfBjCI2Lju0QWPrnc7lxBCCPFVkL6vomrsZtBtSOMBTPrYl4Np1tz5j91EwJNMduhg3l1WFj+loqJqYLrSibZhWfHRbmcQztasWZNnmrHz3M4hnJlmLGhZm/q6nUM4a8+al7avomrCYgLGFBeGniXFCzXTruY/Dy8Pv5Kuvg/XowseexgApk2b1olrKERa3XUw7yfiEVrrHQDk044OTGvcBmCp2zlEy2prM7sDPBnA625nES1jRhGQKAWwye0somXtWfPSdnBBhDvBOAWEU5hpzH4+vmi3TzRuCd2SnchNXDp3/txw07JZU6bkVdvZhUTU4K32rkx2avghMy5l4Ml58+f9euaUmRdyjf4xiO5/9K+PHtSzUJTiMtvm2rZfE9GWiPTf3M4gnBmG2q01lrmdQzgjwluA/sjtHMJZe9a8tP0DHwqFjG3btrXYfns9FXXaddP6KcatIDp97hPzzgaAqVOnejMavO8xYy0RfQLWJ8xd8NhFTe+Zed1N5zLTTd7qzOseCT8iBwlCCCHEQUjbNRfhcNhetWqV3TO35zk9cnrc3Cv76It6ZffK75bZrUd7Pm5daXoMUGc0X5aRyJgAQnzegnmhufPnTgfR0TOumVHQ9HNm4ypAHZfMbVg6c/L0V2deN+2Gg+mzrCx6biRSdVpbrYNID8uK3eZ2BuHsnXc2dTfN6HVu5xDOIpFYUVlZ/BS3cwhn7Vnz0nZapLCwMLMTdXoFzOcQsBOKntbQfTzK+HXxqOLxJctKVqar7+bmLpg3evqk6WeRgYc/W8gYBJD1+Wt6nzx6GIAoAMydP9fpgWsAANOMjQH4DKUwf8iQgjWWFf8es+6dTGb+BKgbopR9oWXF3vT7fQ9VVFQNtG37ekCVBQL5JakL1HgkEZ73+wvejkSi1xLhZMDzaCDQ7yPTjP2YCJl+v++ud97Z1D0zM3E7gGggUPAH04wFAZ5AhFV+f8FLkUj020Q4C1BPBgL5laYZmwFwX683977Kyh61J5wQ+xmz2hEM5v8iEon6iHAjEUy/v+CZsrLouUqhkJleCAZ9b5lm/GpAn2oY9rzBg4/bGInE7gQoNxjMvyMSiXYhwg+JEPf7C35XUbFhiG3Tlcz0RjDoW25Z8dHM+htKeZ4eMqRfRSQSm0bE/evqkvefccZxeywr9gARdvv9BQ+Ulq7v7/GoacxUEQz6nrasqnOY7VFExjK/v/+bkUjsSiIeojX/bujQAXHLiv6IGZ39ft+P3nlnXV5WlucOZtoYDPrmWVbVqcz21UTqX35//tJIJFZExN/SmhcOHTrAMs3ojQB8zHgwGCzYFYnE7we4Ohj0/YwZ400zejSgKgOB/CcjkdjZRHyZ1lgxdGjBPywregUzAkT0B7/fF41E4rOIdLcPP/TdOWjQtuxEovouQH0UCOQ/WlYWP0UpPYkZ/w4GC/7PsqKXMGMEQIsCAV/ENONTAT2gvt77yzPO6LvdsmL3MaM+EPD9tKysqo9S9s3MeD8YLFhgWdEzmXE5QCsDAd/rphkvBvRQwzD+OHhw//WWFbuNmY9OJrf/GOjj9Xjq7yZSm/3+/N+Ul0dP0hqTAXonEPA9Z5obLgbofKVUeMiQ/NWWFfsuMx+ndfLhoUOP32aasZ8Qke33599bWRk9JpHA9wF8EAgUPGGaG04HaCyResXvz3/NsuLjmfVwZv3nYHDgWtOM/QDgnrm5yXs8Hg/t3o05ALYGAgWPmOb6EwA1hYje8/t9z5aVxS9USl/IbDwbDPZ/zzRjUwA+wTD4kcGDB2yNRKJziECBQME9lrWuJ7PxAyJaq3XiNWYaZ5rRE5nVa8Fg/iup5/bw6UT8hN8/4INIJPp9IhxTU+O5t3PnT+1EIucnAG0LBHwPr1698TjDSH6XWa0OBvPDlhU7n5kvJsJzfn/BO5FIdDIRTtLa+O3Qof0/Ns343QB7AwHfj0tLPz7a46m/jVmtDwbz/2hZsaHMXAzg9UCgYGUkEr2cCGfaNi0YNsz3vmnGZwL6WKLsn+7Y0au+a9fYfcy0PRj0/dKyYgXMfAOASCBQsCgSiY4gwiVa0/8NHer7t2XFJzHrU4i8c/3+vptMM3YXM7KDQd+dlZUbuyUSyVkAYoFAwe8jkWiACFcoRf8YMsS3wrKqLmO2zyYynvL7+/8nEolNJ+J+WVn1PzvxxBOrLSt2P4BdgUDBg2VlG/KVopuYqTwY9P2tvHzDt7SmIkC9GAjk/zMSiV1FxIOTSf3Y8OEDqywrdgcz5QUC+T8qLV3f2eNRP2KmqmDQ91h5+cbBWievYqY3AQSU4mGmGZ1kGPy3wYMHlFtW9CZm5CeT+oHhwwfuNs34A0S8x+/33V9Rsa6fbRvTAfWfQCD/KdOMfxPQlyrFy4cMGfCGZUUnMsPPjN8HgwUx04z+EEAXv993xwcffJBbV5d5J0CbAgHfXNOMDwL0NUTGW35//xfKy2OFWvO5zHgmGCwwTTN6A4ACr9fzi0GD+u2IRGI/I0JtIOC7zzQ3HgskZxKp//r9+X8tK4udpRR/mxkvBYMFq0wzOgFAUCn1+JAh+RsikdjtRNx9507fXd26bc1krv0xkfrY78//7erVsZMNg69lxtvBYMHzphkdCeA8Iirx+31lkUj8eiI9MJnMfGj48D7/M83YTwFKBAL5PyktjfX2ePh7zFgTDBbMt6zoGcwYQ0Qv+/2+v0ci8RCRHmbbnj8NG9ZvnWnGbgW4h9dbc/fu3UcZOTnJe5ixJRgs+LVlbTiRma4D6N1AwLckEolfRKQv0NpYrDWPtqzYTmY+nsj+ld9/3CemGb2XGRwMFsw5/H9xP5e2g4s85M1icOck2/0MMooJKAgvC/+luKh4IIEfmzps6kmPr348ka7+HTH3IsI/m14S8WYAxxxME0rx+8mk2pNMehufiaLf1lrlAR8nMjP7/F9DQ2ILwB+luuMdWqtXlcLmxtdRZvWqYRibGl+bzMZHOTnVn6ba1v9gVh4AsG2u0Vq96vFgJwB4vby5vl69yqxijTnW2LZRDTQ0PllWv6u1sQbY1hAK9dCRiHqVCE2ndnZprV41DGq8BwfHtDZe9XoTVBOfigAAIABJREFUGwHAMOzyRMLYUleX3J1qW70BcAYAdOmCup071avAZxfmbtVavUqUjKfWIfmB1v/f3p3HR1Xe+wP/fJ8zM9nYdxBIArhGmJmAVmyr1LYqW68LE9yJeEsrCW1d297WK7b3tr1auwnY0nsValuFuCNQLW3xZ0tdkpkzSYNYhJkERAKyk21mzvP9/TEziiwnQic5k/B9v16+XpnlPOcz55zny+OZZ85xtRK1N6W20Fta8z8HDkQ7EXEw2LBOKbSllj2Q/EyqCQDi8XiDUq51WscbkjmMukRC79I6sT/ZtnpNa+QSEW/evLn90CG1jogOJJuK7dLavU6pRGPysbVZa1fM7bZ2AoDWRjXA7/bvr9uSnwl/AlQs9do3Ae12u61dAOByJbbF4+51gBVNtmX8Q2v+gEjvS+X+m9YqLxCA3rhxd0zrXuuUsg6mcuzW2rPOMKxtAJBIuN4l0omcHP1+ar01iYTaalnckmrrz0rpBADk57ccamnJW0ekU8eSe7vW1joiHUm9t15rtdftpn3J3NjArPInTpxo1dTUQOvB6wxDHwaAWCznA6Xi61wuvT25LV1b4nHWWtP7qbaCWquo1gXNqcd/sSytk38WHNa6dZ1SnLpegfs9rfU6w+CtyWMRGwG1Px53703l+DugCkKhcYnBg0F9+zauA7gZANrbc/a43dY65mQfcLmwNZFQ65SiHan1hrRWjW53XjrHeoBSX6P2ata6fZ1h8N7WVleTx2M9rJQy3G6KpHK8TaQOau3ek9q2r1uW6uXxNCXq6yfqceMa1xElt3NuLva2t3/U9wBEkn3ASM8PMLVW2/PzWw4lP5N+1TDISL32sb4Xj+N9pT7qe4aBTYmEOkwU+wAAiPSblqV65+TsiU+ZMvRjfY+I9lsWHdH3ENVarXO749uTn8kKA8aOWKw1XQP+H7PhBoADB9Cam/tR3/N4sDNZA5J9D0i8o7XRYhjJGkCk39LaeMcwjBgATtWetlSOj/W9RAINgFqndbL/uFy6Lh43mgCd6l/qNa05BwAGDlTt+/Z91PeI4ru0NtYplWgkMjYyG4O11n0Mw2pKbkujmpk3DxyI9tTx8ifDUO3J44MOKqXWpfteIpFoTNYAqyHZtusflqV3E+lUvVF/1Rq5ANgwjFjyMyX7HlF8t2W51xlGYltqW77L7Irn5Oh0DagBeMuBA8l9oRT9mciKA0BbW9shj8ezDrB2A4Dbndgej7vXEelosm2j3rJ4j1JI1QDaoDXlTZkCvXHjnnh7e691hqEPAQCzZ7fWH/U9y3JtIdKW1pyqAVxjWSoCoCXZNv5iWdoCgLy8/EPt7a1H1oD3tLbWGQZHUu/dqLXal5uLvcl14e/MKr++/jxrzJgaPrIGtLV5Pvh436Ot8bhiADuY6S5mtrRWDYaR7nu0HuDucyHLsmmB9YFpgXIACEwL3Dl7euARIHkJ8LJpgZbAjMCFXZWl4paKyZW3Vmz48HF5xTcqy+ffd8Tjx792y3z5uZsQQgiRAZ13nQuglYjPPPrJOMXzAbgUkiNYJxDTWwBdBgC333B7fwBTLFBtptoPhxtmJU9Hi2wWCkV/7XQGYS8YbBxhmhG5/HeWM83o3GAwKpf/znJdWfM6b3DB9CyYvjp7+qzriCkXAK6/4toSsujHAN5rOtzk2CV9m1qaXge4aUF5xZuGR72tGD9d/JvFGbt0rdY8hIjkxmVZj85yOoGwR2R5mFHkdA5hjxnDlWK5cVnW67qa16nXmpg9LfAwA19HchBjITnHo5EVAlUvVZ3UtSM6Q+XcyhH5aDn04GOPHcpku5s3b86JxWJcUlISy2S7IrM2bdrU+5xzzsnovheZxcyqtrYpz+sd1tzxu4VTIpFI7p49e6xJkyY5M49OfCJdWfM6bULnzJkz81esqrorcEXgUXLRRWAMZuItbapt3apVq1o6bqHzLXps0Y6O3yWEcAIzDwBQWlIydDAz1wOoJyLL6VxCiI512uAiz8p5PDA9MFBpflzn45mqqqrT5mJUhw+7b1PKtRfAU05nESfW2pr3EoBLnc4hPo6ZC2LAtEMJfOdAAsP3tuveI/NVU56BvzLzA0T0rtMZxcft3093uVwDqwF0ySUGxKnpyprXiXMu1C8UYwcTLaFm7CybPmtpYEbgtJjwQ0QHmdVhp3MIe0QsZ66yUFsCF0Ra8dD9UZRMCmLwZXWUd2UdRj+9G1fFNB5g5jynM4pj7CMysuKMtDixrqx5nX5/j/Ip5bmtBYenMWM2QDMAbGNgWdWaqh919rqFEN3PwTjf84e9+F75O8hpsZI1yiDwpf1gPXY2agpzcQ8RveZ0TiHEiXXmT1EBAMvWL2tbsfrpZ3e1fHAjK9wEoICAH3b2ep0UDm8fWVe3dajTOYS9YLBxktMZ/hXMTMw8lJkvSP03lJm77IaAnWW/hWGN7UB6YAEAFoP2xEFb2lAAYLiD8cRxhMPR4jfe2D7Q6RzCXlfWvE6bcwEAgUDAow7z51lhFpiugkYuAc+DeXlnrtdpWie+BJDMuchyRPwwuumcC2buGwOubLcw/YN2DCcFHuDGrlwDq5n5D0S0r+NWslNvAwcGe459PpeA4R60A9h/7KvCSVrjhpycuMy5yHJdWfM6bXARmBa4h5rxH0zUF4z/B6Z72txtVS+++GKP/+kfM9drLd8/ZjuleJXTGU4FM6tWC1PrW3H36g9wXs0huBmArzcSMwehxJcHMHMVEXXZDQIzKc+F0IR8xK8ahJwXPgAxgDM8wDVDwIU5+ABAtdMZPylmHthsYaYFlDa1YVB/D5rzDLzjUnghB9jag379EtJabXM6hLDXlTWv006hlk2d9SNANStWTzz1h6einbUeIU43zDzigIXv/Oo93Paz7fC8H0v246Ee8J0jEZs7HI8NcuN/iKjB6ayngpkHHdIo3xfDt6JtyP8gDozNhTUiF+ZgN74N4G9ElPX3QmBmo1Xj7rcO4auv7cfw7e2gvgboor5ovaQ/VuYYuK8XfXivESF6lIyeubh2+rVj3KAznlr99GukjOUMHqahi8umlxUf/d6Vq1f+JZPrziamGfVbFtomTix62+ks4sRCoegNfn/R753OcQqG7I5hyDstUOmBBQA0xUAbW6Deb8eIQW4MAdAtBxcA9vRWeFQZWNc/T1/WZsTG9+2d+0KiBW/Ajd3dYWCRUtRmYc6vd2DUcx9ANVvJO7N96gBc4/Jx3ZgcPA9gtdMhMyEYjE42DOz0eosiHb9bOKUra15GBxcG6DbNVAGgH1h/G8DNNm/v9hPPToQZk10u3gtABhdZjb4CoDsOLiwPQbsUjvlH1k1gj4EEklfE7ZZSg4dmAGYoFNkPwDvEX/y8w7FORcnOOAr+vB/UnJqcygBePwh6uxnGKA986CGDCyJcxszVAGRwkdW6ruZl9Ncih3Tz99pc7aMAIK+lYF6bq73Pif7L5HqzDTOtV8p4y+kcwp5S+LHTGU7Re/1c2D6pF/RZ+R8NMMbmAhf0Boa60QBgu4P5MiYvr30PES1zOscpMjRw7AgQgMUAA8ZxXuqWiLCGyO3Y/aLEJ9OVNa/Hnj0QoqdiZpVI4JLN7bi3rhmXbGqBiwGckw/L2wuvjfPgQZcLr/agyYLdEjOfuy+B5+Zuwrg/7IVq08l6O74AWH4uWs7Jw035LnrO6ZxCdIaMDi6umXrNSLfLPeSTvHfFqhXBTK47m4TDDTO11of8/uL1TmcRJ2aa0Z/4fEV3Op3jVDBzfmsc49sJnzloYRwY6OPBuwUKf3MDdUTUI270VVe3dWgiYczx+wsfdDrLyWJmd0zj/j/vx23r9mLA9hiojwG+pC+sqQPwUr4Hd+YT9YgzTKFQ9AZm9c/S0tHd5pc8p6OurHkZnXPhInUnW/qOT/DWVgD5mVx3NtGaRymFvU7nEPaYaaLTGU4VEbUAeCP1X4+VSFAeoEucznEqiCjOzIsv7Y09pb3gf68dAwe5cbi/B+96FFZ4gJ70S5Fipaw9TocQ9rqy5mV0cNGvacA33x/x/ncBID+RP4FJvwDmB0kZL2utm4jocwB/Dcw9+lQgc7zK5cqVU9JZTilrntMZhD2Pp2VHLJZ/v9M5ThURvc/Mj+QBY4Z4MAjAIQBbUoPDHsOy8FgikegRZ8t6sq6seZ13nYtpZb8GdPvKNU9XHvn87Gmzz2Lof7jIPfT3q3/fba8iKIQQQojj68R7i/AYAh1zNc4Va1ZsBtCquf2szlu3s8Lh6G2mGb3G6RzCnmlGX3I6g7BXXb1ltGk2POp0DmHPNKN3hsPRy5zOIex1Zc3rxHuL0AYGL7huxqyX+7w/4G9La5bGA1MDg4noDoDJalHysyUhhBCiB+q0r0XKp5TnNuc3v0DA5QAOEfAeA2eD0EzMN61Y8/QLnbVuIYQQQjin069zEZgauFQRrmZwARFqNVNV1ZqqnjRL+hivv765T+/efXVJyZDDTmcRJxYMNo4oLR29w+kc4sT+8hd2DRiwZYDXO26X01nEiYVCkX6tra72iy8e1ep0FnFiXVnzOnHORVLV2qpXNcNLRG0rVj/9SE8fWABATo77pkSieYbTOYQ9In7S6QzCXr9+0ZFaGw85nUN0hCry8xOXOJ1C2OvKmteJcy5OX0rRNq25x99avrsj4hqnMwh7Lhe3JhKGzM/KfhGtDbnORZbrcTUvMDXwl9nTA484nUMIIYQQna/TvxYBAGJ6SFO3vPvkKQkGG86rq2sc63QOYS8cbpjpdAZhb9OmTb1NM/o5p3MIe6YZ9YfD20c6nUPY68qa1yVfi6z8w8o1XbGebEHEU7TWewFscTqLODGtcTeAVU7nECfW2pozEOByAH9xOos4MWZMA+LV6CF34+2purLmdfmci9lTZ1/EpJ9cuaaquKvX3VWU4qBlscyaznJE+rQ5m9ZdGYY6oDVWO51DHB8zKwBj9+07ZLlcajgzDwewk4iOd6d54bAeXfOumzFrStm0gBx4QgjRjTHz8GaLf7A/zuG3m3lrpIW37o9zTZvF85i5v9P5hLO6ZM7F6SYYjFwaCjVe4HQOYS8cjt7tdAZh7403tg80zcitTucQH8fMrjaNylf34bab38Z5V4RRNKUWhXdvwYS6FtwXA77odEZxrK6seTK46AREVKKUJRM6s5zWJBM6s5zHE+/NjClO5xDHOKtZY/zK3ej/x30wGttBDa1QT+6CsWE/Bh2M40KnA4pjdWXNk+tcdAKlXC8qFYs7nUPYY6a7nM4g7LW2upry8/WDTucQxxi8vQ29d7SD2vRHV3putkCRNqimGIY7GU4cX1fWvIwOLgLTAlcS8+ft3sOaRmVyndnI6x0pM6a7gdLS0dVOZxD2UpeTlotoZZ/mwR7E8o1jX+jvAvdxQS4imIW6suZl+szFeBCV2b0hNZOzMcPrzSrhcMMsZuugzzfmFaeziBMLhaK/9vuLvux0DnFiwWDjCKWsr/h8xfc7nUV8zDv9DLw3bQCst1tgvNOSPHvx6T7A5L7AQDfqnA4ojtWVNS+jgwtFau2K1Ss6vA/A7KmzbsnkerON1jxEKZKvnLIeneV0AmGPyPIwo8jpHOLjiOhQLMZLZg7EoIm9cemWNuR5iHlcLrUNy8EyI44/OJ1RHE/X1byM/gPIsK4vm16Ws3L1yuPOSA18IdCX3PRLJr4OwG8yue5s0qtX/P9isZj83DbL5eW1ys3lspzPV9RYW9s03+kc4lhuN8y+wPy+LlwwLg/jLc3tbkVv9XKhDsBep/OJY3VlzcvoLdfLppVdC3AVMX6wYm3Vdz/22tSyTxPx7xgYAfD9K9c8/cNMrlsIIYQ4WakLgQ1sTsCrAW+zRp8CwjYDqMl34x0ianE6Y3eU0Z+irlyz8hkCz2PCf5RNm/VdAAgEAsbsqYEHQPwqAy1kqIt6+sAiFIrOD4cj1zmdQ9gLhRpedTqDsBcKRYpCochyp3MIe6FQ9DumGbnC6Ryn4jAwaE8CX97cht++shf//dIHuPfvh/Cz7QkstYAvpAYfPUJX1ryMzwtYsebp/509bRYYtLRsWlk+DvMUJlwE8CP5Lb2+uWz9srZMrzPbENFBZjrsdA5hj4h3OJ1B2FPKnWBO7HY6h+jQPiKjW/4fPiVQur0N936vAb1f2A2lAQzxgL8yHL57RuOu3gaqAfSIWtGVNS+jX4scafa0Wf/OoKUADjBQVrWm6o+dtS4hhBDiVOyP8zdf3ouF121E7pET5cYXwPr9eXj7/AJ8m4hecixgN9Vpp3tWrHn6fwk8D0AfEEZ31nqyUTi8fWRd3dahTucQ9oLBxklOZxD2NmzYlmeaDSVO5xD2wuFo8RtvbB/odI5TEdPotd8CHz0DP8agDxJwA+jtRK7O0JU1L6Nfi1x35azPM9ElHz7BAAj1xFg6e1pgPBgH0i+tWFvVY3+3rnXiSwDtBfCU01nEiRHxwwAudTqHOLG8vMRQZtwLYI7TWcSJaY0bcnLi1QBedjrLycojvFeUCx7qATfFkmfzPQSMzgEXedCCHnRdpq6seRkdXFiKSgkoP85L2xm4+qgvYXrs4IJIv6u1kjkXWY5Ir3c6g7DHnGhWyv2W0zmEPSK8rTXvdDrHqcgB3hyXi023j8D5L30Ao52BM3LA1w3G4RE5MAGEnM6YKV1Z8zptzoUQQgiR7Zg5LwFc2pLA13fFMWxXHO6iXBwa6IKZo/BTIvqn0xmFAACYZtRfUxM91+kcwl4oFL3B6QzCXnX1lr6hUHS60zmEvWAwOjkcjhY7neNUMTMx80Bm/jwzX83ME5nZ43SuTOvKmtdjfr+bTZgx2eVir9M5REfoK04nEPYMQ/UH2PZ+RcJ5RLiMmbvt5fSJiIloDxH9iYieI6IaIoo5nSvzuq7myf0vOgEzrVdKtTudQ9hTCj92OoOwl5fXvqetLXeZ0zmEPSKsIXLL9UiynNQ8IYQQQgjxEdPcenkoFL3Y6RzCnmlGHnA6g7AXDG4ebJrRSqdzCHumGb0mGNwqXwVnua6seTLn4gQW3LigT8UtFaf0HSKzGqcUn1YXDuuOmNUUpzMIe0SuAma+wOkcwh4zzlWKhjmdQ9jryponcy6O42tzbp+oyXqSoP5YeWuFP+aOX7p06dL4J12eOV7lcuVanZlR/OuUsuY5nUHY83hadsRi+T32mjg9hWXhsUQi0ex0DmFPap7DKsorvrHglgX+1N+PV9xSMdnpTEIIIUR3cVp8LXJH4I68yvLKwJHP3Tt3bu+KORVlleWVV90RuCPvyNcWL1v8M6uXtXnBnIolROwd3Dq45mTWFw5HbzPN6DWZyC46j2lG5WZEWa66esto02x41Okcwp5pRu8Mh6OXOZ1D2OvKmtfjvxaZf+v8UXFuvwtEFwKoAoB58+a5W2Lu1wBsBrArnt9eseDWig3MmMHAE4uXLf7ZkCFDWnZFdy1Vmkbtzt99FYCVn3SdWiNHqZ53AZaehpl6zA2JeirDUIpZ5zudQ9hjRh7AbqdzCHtdWfN6/OBCaVoCokE44qZ3nrhnNogbFi9bHACAyvKKECdo2aInFt0PAJVzKh9pamz66aPLHjUry+evBdFJ3e3P7y9aktlPITqD318oNy3Lcn5/cRRy07Ks5/cX/bfTGUTHurLm9fjBxaLli2dW3FIxmQw8/OGTjBKAwh89prfJpScCiKSeecrQtGxBeUWUGf3jnvhNR7dbXb1ltMfjGhyLJd6dNGnsgZqa6LluN+Xv2TM63K9ftJdSxngA+7ze0f/YtGlT71gs7yzL0nv8/uJofX1kmGWpMxKJWGNp6Zm76+oaxzJzv1jM/fakSSNagsGtXsNQhtdbFNy8eXNOa6vnfK0Th3y+sf+sr982wLJ0sdZqp8836r26undHMbuHWJbe4vcX70/niMV2106cODFRV9dYykztXu/of9TX7+plWa1nM/Ner7coks7BHN/m9Y7blc7BnLPJ6x3WXFe3dYLWyuX1FgXr6+s9ltVrPLN12Osd805tbUN/AGOYjSavd+T2dA4AWydMKNxXWxs5B1AFsdju2kmTJsVraxsmAio2YcKounB4ZwFR+zkA9k2YULi1rm7rUGZjJJG1ffz4MU21tQ1jAPQ3jLx3SkqGHK6t3TYe0J4JEwprqqur3R7P4AmAbp4woXhTOgdRfNf48eO2hcPbRxJZQw1DRUpKRu0Nh7eeTWT0MozDdSUlJbFwOFqqlE6MHz+mtqZm2zi3W/clov3jx4/eEg6/O4TIPcow9HslJcU7w+FoMRENSOcIhxvPJ+Kc8eNHB2tqalwez+AJ8Ti3TJxY9HYoFOlnGGpsOodpbjtDKT0sncM0t5yllKt3Xl7sH2eeeWZ7OBwttSxtlZaOCVdX78j3eOLnxuPqwMSJo94NBjcPdrk8o9M5QqFIkWGogR5P6z/POeecQ0fkCK1fDzVwYKM3kUBraWnhxurqLX09Hte4WCyxe9KksY3BYOMIl4uHt7Ya0U99auSeUGjLmYbh6tO7t64vLi5uM82on5nZ7y82N2zYlterlz4vnaO6escgjydeGIvxjkmTit4PBrcWulzGoJaW2OaLLjrzoGk2lCiF3E2bRpsAcM45jb50jtdf39wnP99zZiJhfVBaOqYhnSMWczdMmjTig3SOw4fVxosvHtUaCkV8AOD3F5uRSCT30CFVYlmJg/v3j4kMHBgdB6iCRILeLy0dvSOdI933gsGG81wu5G3aNNoMBMB1dY1+rdHm8xXWH933qqujwz0eGpHue+lj4Mi+p5RSPl9R6Oi+98Yb2wfm5VlF6b6XrgFH9709e0aHp0yBdWTfS+foqO+lcxzZ99I50n2voxpwZN9L14B030vXgKP7XjpHuu8dWQOYDfeECYU16RpwdN9jNpq0jh/W2jXI7dZ90zXg6L5XW9swkciKjx8/pvboGpDue+kacHTfS9eAI/teOkdHfS+d48i+l64B6b7XUQ04Tt8Lrl8PY+DARm+6Bhzd99I50n3v6BpgmlG/1lqXlo4JH9330jUg3ffSNeDovjd+/OhQVRXoyL6XznG8vqdUwkMULziyBhzZ9zL5b+9pMefiGMxDCRxNPyTi9wF8+DOqRcsX/W3RsiWXaI07Fi1fPHPp0qUHjm7C5TImaq0DHo8annxMl2utA717v+9h5psB625AXwkA8Xj+EK11gEhdmHxM52qtA4bhGgsAiYT1Ga11wOWy+gOAUupLAK4BgL1783pprQOA+lxyWasw2VZifHJZozTV9hkAYBj4vNY6kJc3OqeqCkprHdCapwFALNY8KNkWfwoA2ttxdvJ11zgAsCyerLUOKHV4QLJtNVNrmpVcNj8/uSxdltxm1qjksokJAKC125fKOTLZFl2mtQ4kEn3zmJm01gHmxIzkZ2geqLUOWBZPTi7rOiuV46xUjou01oH29sMDk7srMSO5PFMi0TcvuSx9PvkanZFc1u1LtpWYkFw2Pir1+ue01oFYLD8/+TrNSiTUzOS2sn6bynFxKse45LI4O3WgfCr5GVoGJ1/naVrrQFUVVF7e6JzkPqQvJNejRiTfa5Qmt09ivNY6EI9bhakjZorWOrB3b16v1CF0TWo/Q6lEv+R2T3w6mcs1NrksnZtsS12YfJw/JLmovlJrHaipqTF6937fk8yBLwKAx6OGJ48lY2Kq7fOTx4NVlHxsXKq1Duzbl0ifHr1aKfVvANC3b6Jvsq3EZ5PHdNsYrXXA7abzkrmMC7TWgZwcNTS1fa7QWgf8/nddF1wQdSc/g74CAHJz1bBkW2pS6r3npdouTn4m47PJ/hLvk8pxFRGuBoCDB60+ybaMS/r1i460LPpx8hjnklT/mJjKNSzZlr5cax0YOXK7Z+PGja7UMX4lADQ3FwxNHZcXAIDbjXTfG5PaPp9Jbft+qe3zb+m+d+BAQe/kcee6NPmZEoXJthLnJ7eHKk0eezQ8+Zi+oLUODBjQlLN+PYzka3pqsq/lDkodwxcm+xPOSR57yRqgNX86+Zn0gGQfUF/SGtcm35tbkFyWPpd8nBidzGyl+p7Hn9pe6b73ea11oFevXrkAKNVPpydzHE73vYuSx22yBjC7031vcvJYax6YbFvN0BqzAKC1NTfd91JzK9TIVNtegCoMw7o+uZ+sUcnXkzWgtTU3L9kWZmmtZiS38+EByWVxcTKH+6xkjg/73kXJfts8KLlsYnpy24N69eqVm9ruX0j1j1QN8PiTbVsTkssmRifbpinJx7kFqRzXWlay77lcVv9Ujfx08vMn+14shnOSy/KFycd5qRqgp2qtA+vXwxgwoCkn2dfoi6ljekTqmC9NfoZk38vNTRQm23JdqrUOHDhQkO5716SON3g8VroGfCZ5LLlSfQ+pe1SpC5KPc1M1gK/UWgc2btzoGjlyuye53fXlyWOchiXbUhOTubgk1XYRET+plHGJ1jpw8KDVJ/k6rgZwFcTJq7ilYnLlrRUbPnxcXvGNyvL59x3x+PGv3TL/c5laXzjcMDMUikzJVHuic5hm9CdOZxD26uq2Dg2FGu51OoewFwpFbwgGGyc5nUPY68qa1+O/FjkeYnoLhP8C8P3bb7i9P4ApFujuTLXv9RauylRbovP4fEV3Op1B2Bs/fkwTgAedziHs+f1Fv3c6g+hYV9a80/JrkaaWptcBblpQXvGm4VFvK8ZPF/9m8Z5MtV9X1zi2ru7dUZlqT3QOObuU/cLhnQWmufVCp3MIezU10XPr6yNyhc4sJzWvi1TOrRxx79y5Gf9pTigUnR8OR67LdLsis0KhhledziDshUKRolAostzpHMJeKBT9jmlGrnA6h7DXlTXvtPxaJG3RY4t2dEa7SnHQsri1M9oWmUOk5VRuljMMdUBrrHY6h7BHhA2Afs/pHMKe1DwhhBBCCPGRYDByaSjUKHdyzHLhcDRjk3hF53jjje0DTTNyq9M5hL1QKDotGGw4z+kcwl5X1rzTckIW3ggdAAAdQUlEQVRnZyOiEqWssU7nEPa0pplOZxD2PJ54b2ZMcTqH6JBfKS2T2LNcV9a803rORWdxudTLQCzmdA7RES238s5ySuXtZm5Z7HQOYU9rPJubi31O5xAdkZonhBBCCCHSwuGGWaa59XKncwh7oVD0105nEPaCwcYRphl5wOkcwp5pRucGg9HJTucQ9rqy5smci06gNQ8hogFO5xAdobOcTiDsEVkeZhQ5nUPYY8ZwpbhPx+8UzpKa161t3rw5p76+3uN0DmFv06ZNGb+AmsgsZlbh8M4Cp3MIe5FIJLe6utrtdA5hT2qeEEIIIYT4iFz+u3uQy39nP7n8d/cgl//uHrqy5smci05ARAeZ1WGncwh7RNwpl38XmaOUO0FEu53OITq0j8hocTqEsCc1TwghhBBCfCQc3j6yrm7rUKdzCHvBYOMkpzMIexs2bMszzYYSp3MIe+FwtPiNN7YPdDqHsNeVNU++FukEWie+pDV9zukcwh4RP+x0BmEvLy8xlFnf63QOYU9r3JCTE5fBepbryponl//uBET6Xa1lzkW2I9Lrnc4g7DEnmpVyv+V0DmGPCG9rzTudziHsSc0TQgghhBAfMc2ov6Ymeq7TOYS9UCh6g9MZhL3q6i19Q6HodKdzCHvBYHRyOBwtdjqHsNeVNU/mXHQCZkx2udjrdA7REfqK0wmEPcNQ/QEuczqHsEeEy5hZLi2d9bqu5smci05AhL8nEtTmdA7REf6V0wmEPcvS+wzDWOl0DmGPGX9WimTORdaTmieEEEIIIdJMc+vloVD0YqdzCHtyK+/sFwxuHmya0Uqncwh7phm9JhjcKl8FZ7murHky56ITMKtxSvFop3MIe8xqitMZhD0iVwEzX+B0DmGPGecqRcOcziHsdWXNkzkXnYA5XuVy5VpO5xD2lLLmOZ1B2PN4WnbEYvn3O51D2LMsPJZIJJqdziHsSc0TQgghhBAfMc2Gm0yzYYbTOYQ904w+6XQGYS8c3j7SNKM/djqHsBcKReeHw42fdTqHsNeVNU++FukEzNxHKU44nUPYY6YRTmcQ9rSOuwAMdjqH6FB/Zivf6RDCntQ8IYQQQgjxkddf39ynvn5XL6dzCHvBYKOM4rPcX/7CrnD43SFO5xD2QqFIvw0btuU5nUPY68qaJz9F7QQ5Oe6bEolmmXOR5YhY5lxkuX79oiO1Nh5yOofoCFXk5ycucTqFsNeVNU8GF51AKdrFzHudziE6wv90OoGwx2zEiBB1OoewR4T3taaDTucQHZGaJ4QQQggh0urqGsfW1b07yukcwl4oFJnidAZhLxzeWWCaWy90OoewV1MTPbe+PiJX6MxyXVnz5GuRTpBI6Cu0Nj7tdA7RESX3FslyWrcOZqYKp3MIe0rhmngccm+RrNd1NU+uc9EJlOKgZXGr0zmEPSL9e6czCHuGoQ5ojdVO5xD2iLAB0O85nUPYk5onhBBCCCE+Eg5HLpLbD2c/04x8xekMwl5tbUN/02woczqHsGea0c+Z5paznM4h7HVlzZM5F51Aayo1DDrX6RzCHrO6wekMwp5l6b7MerrTOYQ9ZlwMqGKncwh7XVnzZM5FJ3C51MtALOZ0DtERLbfyznJK5e1mblnsdA5hT2s8m5uLfU7nEB2RmieEEEIIIdLC4YZZprn1cqdzCHuhUPTXTmcQ9oLBxhGmGZGfDGc504zODQajk53OIex1Zc2TORedQGseQkQDnM4hOkIyAS3LEVkeZhQ5nUPYY8ZwpbiP0zlER7qu5smci07Q3h7/be/efbXTOYQ9Zrre6QzC3v79RdsHDNhyj9M5REd4cUuLq93pFMKe1DwhhBBCCPGRUCg6PxyOXOd0DmEvFGp41ekMwl4oFCkKhSLLnc4h7IVC0e+YZuQKp3MIe11Z82TORScgooPM6rDTOYQ9It7hdAZhTyl3goh2O51DdGgfkdHidAhhrytrHnXVioQQQgjR9ZjZADABwAgAuwCYRBR3NtXpjSrmVJSc7ELB4ObB9fXb5NciWS4c3nq20xmEvfr6ek8oFClyOoewV10dHf7665vl1yJZiJnHHk7wk3tjXB06qDftjnF1s+bnmPn8zlyv/FrERuWtlXdB6ysAfPFkliNyBxKJxF4AT3VOMpEJWhtLAVzqdA5xYrFY/ggADwCY43QWcWKGgbkul6sawMtOZxEfYWZPG+PRZz7AxUt3IO+QRZSvoG8aipLbRmAIM08looOdsW4ZXJzA1+bcPpE1n82nsCyRfldrmXOR7Yj0eqczCHvMiWal3G85nUPYI8LbWvNOp3OIY1xwKIELvx9F/pY2gBkEwNjeDjV1ICaMycUlAF7qjBXLnIvjuPvmuwvaXK1PaOb5ivHEouVLTurMhRBCCOE0Zi5/8yB+cXEIvSz++L/3q8aj9cr+uN9t0EOdse7T4szFHYE78uIF8RmLli2qSj9379y5vZutvKlEFHM3u19O9Ip9ixkzGHiinVvPIsZLBvHZDNWn8ubK4kVPLIp80vWZZtRvWWibOLHo7c75RCITQqHoDX5/0e+dziFOrLp6S1/DMD7j9xetdjqLOLFgMDrZMLDT6y36xHVSdInDAzzHP4sw2A2tFDrlKxHgNBhczL91/qg4t98FogsBVAHAvHnz3C0x92sANgPYFc9vr1j0+JIvArgfACpurZijGGcAqhDAQHbxBACfuNMwY7LLxXsByOAiq9FXAMjgIosZhuoPcBkAGVxkMSJcxszVOIk6KbrEhmFutN80DAW/2Qno1NmLqwcBZ+YBCljfWSvu8YMLpWkJiAYB/OH0CU/cMxvEDYuXLQ4AQGV5RejIsxOLH1+8HADmz5/fSzXj4sXLlrxwMuskwt8TCWrL5OcQnYF/5XQCYc+y9D7DMFY6nUPYY8aflSKZc5F9duYZuOfekXhwRn/03toGY3QuEiUFiPV24z4AWztrxT1+cLFo+eKZFbdUTCYDD3/4JKMEoPBHj+ltcumJOGrUvWTJksM4wS9FwuGGrwM8A+D7vN7i102z4ZdEPLa1NX5tQUF8Y0uL+6VwuOE9r7ewPBjc6jUM9WOt6Y9+f+GD4XDkeoDmAvQLr7dwlWlGv0eEyVpb8/3+sZtNM/o0M/L9/qJpwWDjCMPQy5lR6/MV3RUON3we4G8BeMrrLfo/04xWEuHftMb9fn/RhnA4uhjAWURUtmfP6EMDBjSsZcZOn6/o5trabeOZrZ8A/Gevt/iHptlQRsRf1poX+/3Fz4fDDfcD/BlAV3q9Y94xzcgKQPXz+QqvqKvbOlRr9Vutud7vL/5GKBSZohR9h5mqfL7CpeFw5HaArgHU97ze0a+ZZvQXRDg3HvdcP3Hi8D21tQ2vaI3dfn/RDcFgw3mGwT9nxv/z+Yq+Hw5HrwXwVQC/9HqLnjHN6HeJcKll0ddLSws3hkLR3yuFwRMmFF5eU/P+QLc79iQzNvl8RQvC4cbPAvo/AX7W6y1+1DQb5hFxQGv+b7+/eH043PBTgM9XSt80fvyYJtNseBnQ+32+4tlK6ZpwOPpHrflvfn/xwlAocpVSVMFMv/b5CleGw5FvA3QZkXHnhAmj6kwz+gQRhu3dWzh14MDG3sy8Umts9vuL5odC0YuVwgPMeMHnK1oUDkdvA3AdQD/yegv/ZJrRHxPBa1lqTmnp6B2hUHQNEVp8vqJZNTXbxrlc1qPMeN3nK7ovHG6YCfDXAH7M6y1+MhRquFcp/qLWfI/fX2yGww3LAD4jPz8+Y88e5OTluZ9h5q0+X/FXwuHIpwD6L4Be8noLfx4KRcqVohuZ9UM+35hXTDPyP0RUCrhu9XpHbjfN6CoACZ+v6OpwOFoMYKnW9JbfX/gfoVB0ulL4BjMt9/kKf2ua0buIcCURfWvChMIa04z+LxEKm5uNLxUUuBXQ9jwzGny+on8PBhsnGYb+ITPW+nxFPwmHG24B+GYi/GTChKK1oVDDD5TiC7TmL/v9xdFwuOF5Zq18vuIvVVdvGe12G//HzEGfr/ib4XDESu2n3/r9xctDocg3lKLpzPo7Pt+YN00z8isiGuNy5V+dn9+cOHiQVjHTdp+v8FbTjPqJ8KDWeMXvL3ooFIreqBTKmennPl/hS6YZ/T4RLlJKfXX8+NFbTDP6DIAcn69ohmluO4PIWsaMsM9XdHco1PBFpfheZjzp8xU9Fg43fA3gmZaF/ywtLfp7KBRdohTOjMetWcC+Frd70BqteYffXzynrm7rBK3Vw8z8J5+v+EfJK/jSbVrzIr+/+AXTjDxARBczWxU+39h/hsPRlczo7fMVTa2vjwxLJOgJZtT5fEV3hsPRywB8W2ta6fcX/joUilYohauY6QGfr/CvphldRISzXS5j9nnnjdxfW9vwMjOafL6im0yzoYSIf6Y11vv9Rf8dCjUElOJ5zHjU5yt61jQb/pOIP5tI4GsTJxa9bZrRJwEe6PMVXx4Mbh5sGO7fM/NGn6/467W1Wy9hVvcx8zM+X9EvTTPy1XA4eq1l8X+Vlha/apqRnxPReZYVv6G09Mzdphl5BaA9Pl/R9bW1kXOY6RFmes3nK/yeaUavIcLtWtNSv7+wKhSKfkcpTAHUHV7v6H+YZvS3RBg6YULhFRs3bu+XSFgrmPGOz1dUGQo1fFopXqg1nvf7ixaHQg1fVorLAPzQ6y36s2lGf0KE8S4X31xSUrzTNKNriXDI6y0qC4W2nKmUsYSZN/h8xfeHQpF/U4oqAf4/r7f4KdOMfIuIPm9Z+u7S0jHhUCiyXCkaEY9/MA3on+92G08z07s+X+Ht4XDkIoC+D9Aqr7fwF6YZnUuE67WmB/3+wj+GQtGHlIKP2Sj3+Ua9Z5rRlwC0+3xF19bWNoxh5l8BeMPrLfquaTbMIOKva41lfn/R70Kh6D1K4XIA3/R6i4Km2fA4EY/s04dntrQUuBKJlucARLzeonmhUOMFSukfaM2riehntbUNOfkFOQ2fGdi/9oK+nsPGoWZv+wctvbY16fVjxw6Lm2bkRSKlvd7CqzLzr+5ppOKWismVt1ZsSD+unFPx2ILyirnpxwvKKx5ecOv8ykytzzS3Xh4KRS/OVHuic8itvLNfMLh5sGlGM9Y3Recwzeg1weBWr9M5hL2urHmn5eW/mVDL4DPSjzUwgCzUZ6x9VuOU4tGZak90DmY1xekMwh6Rq4CZL3A6h7DHjHOVomFO5xD2urLmnZaDC2J6C6DLAOD2G27vD2CKBarNVPvM8SqXy/VKptoTnUMpa57TGYQ9j6dlB1ITrUX2siw81taW+LvTOYS9rqx5PX7OxfE0tTS9PrRgcNOC8oo3GRhNjB888pvFezLVfmnpmXKjpW7A6x3zjtMZhL2SkpIYgKjTOYS9SZOK3nc6g+iY1LwuUjm3csS9c+f2znS7ptlwk2k2zMh0uyKzkpPVRDYLh7ePNM3oj53OIeyFQtH5ycnVIpt1Zc07Lc9cpC16bFGn3H6WmfsoxYnOaFtkDjONcDqDsKd13AVgsNM5RIf6M1v5TocQ9qTmCSGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgiRaYFAwKi4peKsQCBgOJ1F2PvGTfOGL1iwIMfpHOLE5s2bl18xp6IEx7+btMgCC25c0KfiloqznM4hPhFK9aeMOS2v0NnV7rjtjgFDC4b8gxQeGFow5M358+f3cjqTOL6vz7397ITLvY4OkVy+PUstmLNgrCfmriPiysry+a8umLNgrNOZxMd9bc7tE9ltVZOir1feWrFh3rx5bqcziROrvLXyLgL/LJNtyuCiC8QT8dsJ/JNFyxZfT4SXjGa+1OlM4lhfvfmrQxJa3cGA/F4/q1lXgeneRcuW3E7AE0xWwOlE4uMsUp8lbcxetGxRBTPecbe5JzmdSRzf1+bcPpE0n53pdmVw0RWIz9OUvHeJ1qhnolKnI4lj/fKJX+5avGzxVwlsOp1FnNgjy5c8vGj5omfmz5/fi0EB0safnc4kPm7xssU/s3pZmxfMqVhCxN7BrYNrnM4kjnX3zXcXaKW+YxHfl+m2T+srdP6r7gjckRcviM9YtGxRVfq5e+fO7d1s5U0lopi72f3yT6t+2spATJFqAwAw4lCq1bHQp6FPup+czChObj8tmLvgC9yi7wPRDx75zSNvOpf69HIy+2jIkCEtu6K7lipNo3bn774KwErHgp9mPul+aletDxHjJYP4bIbqU3lzZfGiJxZFMpFBzlycovm3zh8VL2j/IYjvSD83b948d4vOew1AAMAX4/ntLwIACG9qrSek/h5PWmfsDqzC3kntJ+GYk9lPC269/VLW1r+7mz1XLnp80ctOZT7dnMw+qpxT+UhTY1PRkmVLTBCvhcJAp3Kfbk5mP2mFNwAMA9RlAAayiydkKoecuThFStMSEA0CmNPPeeKe2SBuWLxscQAAKssrQpU3VxZr1s8pomcr5lR8iRQXDCwc8gPnkp9eTmY/ZWrELk7eyewnZtwIYFyiILZqQXkFQFz1yONLfuVY+NPEyewjAE8ZmpYtKK+IMqN/3BO/ybHgp5mTqnmPL1oOAPPnz++lmnHx4mVLXshUDhlcnKJFyxfPrLilYjIZePjDJxklAIU/ekxvk0tPXPL4kqcXLlz4mf3v7hj6s+VL5dbEXehk9hOACAAsWrbk6i4Pepo7mf30yONL5jmR8XR3kvvoaQCXVNxSMXDxbxbvcSDuaetUat6SJUsOA/hiJnPI4CKTmIcS4a/ph0T8PoBhALBw4UINQAYW2cBmP4ksIvsp+3Wwj2RgkSUc6Esy5yKDmFDL4DPSjzUwgCzUO5lJHEv2U/cg+yn7yT7qHpzYTzK4yCBiegugywDg9htu7w9gigWSyZtZRvZT9yD7KfvJPuoenNhPMrjIoKaWptcBblpQXvGm4VFvK8ZP5bRg9pH91D3Ifsp+so+6B9lPPUTl3MoR986d29vpHMKe7KfuQfZT9pN91D3IfhJCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgiHkdMBhBDZ6erLrx7idhnfcnPOwt+t/d1BJ7MEpgYGK8I9TKRWrl55t5NZOrJw4UK18a1/fJ8T9HjVy1XvOp1HCCcopwMIIbKTS7n6A3RHq24tcDqLUvhPBuYwc5PTWTqyfv16Bab/gIFip7MI4RQZXAghsh4zzibCM1Vrqh5yOosQomMupwMIIT6ZQCCQRy34DUH9J7P1ZQZdTMB2Bh6oWlNVd+PUG/vEVez/LKXve2bVM5sA4Nrp144xoP6H2/Hv6I8YteA3THhIacxi4GKA/sYt/H1VwLcy0ywAh5jxUNXaqlfT63UZrjPLppX9nMFnKyAcZ+s/nl377Pb062VTy74M4msA5IHpT4dw6MG1a9e2pzIb1IKnmPHfijELCn1XrK5acJzP5qFmug/gKwC4QfSnNqPtgRdffPFQ2fRZ3wPjfDBGzZ4W+OGKNVXfPnr5a2dee45hGbeA+QoQdhF48Yo1T78EAIErAgNI8b+D6IsAmAl/tbT12LNrn92e3qYns02uueKa4W7D+B+l9Y8tZXwF4E8BVGu4jG8/+eKTxz2zYreN7LIL0V3JmQshuon29nYXGLOY9ctgRQR6AsCnCFgNADHEcsCYpeJqUHoZl3L1A2OWK9+Vm16eNKrA6MuEVwD+KuWjjpluJuKnwXCRwktTp07N+WjN/DSg9wBYzsAXXGT8rXxKeS4AzJ4WeBjEPwJxkIFlTDyrF/X6I1LzuXbv3k1gzCLgd0yYBebw8T4bNeNZgK9lxsNE/D9g/nxuImfdwoULFbR6C8A+BqJa4c3jLW5Y6hkAfQFaCPA/GfTMdVOv9SdfxG9Y0S0EXkHMy4lxtYuMZ47cpiezTQwy+jBws6XUajDeA9EjAE+wrMRfA4GA5+hwHWwj2+xCdFdy5kKIboeWr1y78j4ACFwZCJPCa4FpgWFgWJ9oacayFWur7geA2dMC4xi4kYHJK1c/vTNwZeAVYmzqi/zzLKAluQDds3J11XIACMwMPEcWtjbnN18XuCLwVwa+TsQ3rVj99FMAcM3Ua9a5yNhWNr1sysrVK//y4UoZ76xcW3UtAD46T2Bq4FIA05kwqWpNVU2qnb+6yNhWX13/b1Vrq54rmxb4GgF1K1+qeu7o5a+bdu2FGjiPwTdWra0yAawKTA1s0mT0BgAC9jJw/4o1T9ek1pcgwm9PZZsACH30kbC4au3KHwHANVdc84rLMKJoxmwAT3742a4IjLPbRoqtFrvsQnRXMrgQopsh1q+l/za00aCVBWWoETqht32yFvj1D/8CfwBQqGpN1c7UE7sBQCtjICxuAQBlqD+m31+1qipSNi0QIeZz4FIJMCvWakxgWuDOI1awj5n9AD4cXJDiJ3GcgUXyRVwI4EDJBSWhqtVVAIBn1z67vWxaoBGM8QCOGVAcqc1jveOJqQgx1pVNn/UsWK3jGP9+5bqVBwBgxZqqW2bPmO0tm152N1gXAZiJY87afsJtcgTFan3672dffvb9smmBdxRw3sc+mosusttGbR7rMbvsQnRX8rWIEN2MZmpJ/x0zYsf/Bzv93oTud5zlDx31VLNdGwT62BkRBg6BFIExGIBFiocrQnH6PyL8TjG/fdQ699u0PwhA48KFC/XHnwcrDcMuGwA8//zz+7kFE4j4HmYaAvD/kgfvzZ4x2wsAZdMCL7HWL0Pr8xnUAKIlR7dxstsEAIhp59HPMbH++BP226ij7EJ0V3LmQoieIo4YPAAZlHvEs2f/q81alnUegCYAuOqqq/pRDOeD8XNW/AExXCqhf/nky8/Up95OgemBuQrGP09iFREAJYErAgOqXq7aCyQnTTJQyIR3Olr4+iuvnZBwYeyKl55+HMDj5VPKc1vym1+E1vdeN+3aX2hgmmI98am1z4SA9OTKk9oEx2Up60wAUSB5TRAA55PGj498DxNvsdtGdtkB3PivpxTCGTK4EKKHqFpXdaBsWmAbmL969eVX16ocNYAs3PuvtsvMP7z68qu/5Orr2kMtfCeA/W3utipmtvKs3G2WMr4bmBG4283ulgQn7gF4QVzHz/qk7Rsu4znLSjykFL4O4H4A5HK57gLzrjaj7fmOA8Iijd8FpgamVq2t+n/Ng5s1NaMvM70GlysOSxPDlQMAgamB80B8NwBauHChCgaDp7hVAALuD3wh8Gaid8Lljhn/BcY+TfTKke9pU23r7LZRDij3hNmF6MbkaxEhehAGvgPgS26Xa4dhqZcYeORfbZNAf3a7XI3UjIPQNF8pvvbFF188tGrVqhZWmAXii0gjkuD4ewSeDdCsZ19+9v1P2v6TLz7ZRMzXM2FB2bTAnrJpgX3QfD0rXLVq1aqWDpd/+Zl6AlYQYX3Z1MB2asZ+AL0Ml/HoilUrgmBezqRfLZsWaCLCS2C6G8CujW/U//Vf2y7YRB7sdLe7thHTFCbc8uE8jZSOtpFd9n8lmxBOk8t/C9HDlE8pz23Naz1jxdoVWzLVZmBKoBfyUby7Zffb69evTxz52pQpU1xDew0tQQKevrv7mktrlsZPZR3lU8pz2/IPnc/katb5+p9VVVWf6NcvH2acGSgmTUWGZe1J9Fb1Ry5/9eVXD3GRq0/6ctzlU8pzD+Yc7H8yg6AP13Nl4GxS2GSRHhs34rsL2t2jj/jK47g62kZ22YXojmRwIYQQJ+HIwcUzq5/Z6nQeIbKRfC0ihBBCiIySwYUQQpyEkotKNrvZ03f8BeOjTmcRQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQgghhBBCCCGEEEIIIYQQQvyL/j/QXF0PrcGUKAAAAABJRU5ErkJggg==", "image/svg+xml": [ "\n", "\n", "\n", " \n", " number of samples\n", " \n", " \n", " 100\n", " 101\n", " 102\n", " 103\n", " 104\n", " \n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", " 10-4\n", " 10-3\n", " 10-2\n", " 10-1\n", " 100\n", " \n", " \n", " KL-divergence\n", " \n", " \n", " KL-divergence between exact and moment matched posterior\n", " \n", "\n", "\n", " \n", " \n", "\n", " \n", " \n", " \n", "\n", "\n" ], "text/html": [ "\n", "\n", "\n", " \n", " number of samples\n", " \n", " \n", " 10-5\n", " 10-4\n", " 10-3\n", " 10-2\n", " 10-1\n", " 100\n", " 101\n", " 102\n", " 103\n", " 104\n", " 105\n", " 106\n", " 107\n", " 108\n", " 109\n", " 10-4.0\n", " 10-3.8\n", " 10-3.6\n", " 10-3.4\n", " 10-3.2\n", " 10-3.0\n", " 10-2.8\n", " 10-2.6\n", " 10-2.4\n", " 10-2.2\n", " 10-2.0\n", " 10-1.8\n", " 10-1.6\n", " 10-1.4\n", " 10-1.2\n", " 10-1.0\n", " 10-0.8\n", " 10-0.6\n", " 10-0.4\n", " 10-0.2\n", " 100.0\n", " 100.2\n", " 100.4\n", " 100.6\n", " 100.8\n", " 101.0\n", " 101.2\n", " 101.4\n", " 101.6\n", " 101.8\n", " 102.0\n", " 102.2\n", " 102.4\n", " 102.6\n", " 102.8\n", " 103.0\n", " 103.2\n", " 103.4\n", " 103.6\n", " 103.8\n", " 104.0\n", " 104.2\n", " 104.4\n", " 104.6\n", " 104.8\n", " 105.0\n", " 105.2\n", " 105.4\n", " 105.6\n", " 105.8\n", " 106.0\n", " 106.2\n", " 106.4\n", " 106.6\n", " 106.8\n", " 107.0\n", " 107.2\n", " 107.4\n", " 107.6\n", " 107.8\n", " 108.0\n", " 10-5\n", " 100\n", " 105\n", " 1010\n", " 10-4.0\n", " 10-3.5\n", " 10-3.0\n", " 10-2.5\n", " 10-2.0\n", " 10-1.5\n", " 10-1.0\n", " 10-0.5\n", " 100.0\n", " 100.5\n", " 101.0\n", " 101.5\n", " 102.0\n", " 102.5\n", " 103.0\n", " 103.5\n", " 104.0\n", " 104.5\n", " 105.0\n", " 105.5\n", " 106.0\n", " 106.5\n", " 107.0\n", " 107.5\n", " 108.0\n", " \n", "\n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", " \n", "\n", " \n", " 10-9\n", " 10-8\n", " 10-7\n", " 10-6\n", " 10-5\n", " 10-4\n", " 10-3\n", " 10-2\n", " 10-1\n", " 100\n", " 101\n", " 102\n", " 103\n", " 104\n", " 105\n", " 10-8.0\n", " 10-7.8\n", " 10-7.6\n", " 10-7.4\n", " 10-7.2\n", " 10-7.0\n", " 10-6.8\n", " 10-6.6\n", " 10-6.4\n", " 10-6.2\n", " 10-6.0\n", " 10-5.8\n", " 10-5.6\n", " 10-5.4\n", " 10-5.2\n", " 10-5.0\n", " 10-4.8\n", " 10-4.6\n", " 10-4.4\n", " 10-4.2\n", " 10-4.0\n", " 10-3.8\n", " 10-3.6\n", " 10-3.4\n", " 10-3.2\n", " 10-3.0\n", " 10-2.8\n", " 10-2.6\n", " 10-2.4\n", " 10-2.2\n", " 10-2.0\n", " 10-1.8\n", " 10-1.6\n", " 10-1.4\n", " 10-1.2\n", " 10-1.0\n", " 10-0.8\n", " 10-0.6\n", " 10-0.4\n", " 10-0.2\n", " 100.0\n", " 100.2\n", " 100.4\n", " 100.6\n", " 100.8\n", " 101.0\n", " 101.2\n", " 101.4\n", " 101.6\n", " 101.8\n", " 102.0\n", " 102.2\n", " 102.4\n", " 102.6\n", " 102.8\n", " 103.0\n", " 103.2\n", " 103.4\n", " 103.6\n", " 103.8\n", " 104.0\n", " 10-10\n", " 10-5\n", " 100\n", " 105\n", " 10-8.0\n", " 10-7.5\n", " 10-7.0\n", " 10-6.5\n", " 10-6.0\n", " 10-5.5\n", " 10-5.0\n", " 10-4.5\n", " 10-4.0\n", " 10-3.5\n", " 10-3.0\n", " 10-2.5\n", " 10-2.0\n", " 10-1.5\n", " 10-1.0\n", " 10-0.5\n", " 100.0\n", " 100.5\n", " 101.0\n", " 101.5\n", " 102.0\n", " 102.5\n", " 103.0\n", " 103.5\n", " 104.0\n", " \n", " \n", " KL-divergence\n", " \n", " \n", " KL-divergence between exact and moment matched posterior\n", " \n", "\n", "\n", " \n", " \n", "\n", " \n", " \n", " \n", "\n", "\n", "\n", "\n" ], "text/plain": [ "Plot(...)" ] }, "execution_count": 3, "metadata": {}, "output_type": "execute_result" } ], "source": [ "using Turing, Distributions\n", "\n", "#= The parameter values of Θ_1, Θ_2, Θ_3 follow [Aït-Sahalia, 1999],\n", " [Sun et al., 2013] and [Del Moral et al., 2015]. =#\n", "const Θ_1 = 0.0187\n", "const Θ_2 = 0.2610\n", "const Θ_3 = 0.0224\n", "\n", "#= [Del Moral, 2015] observe the OU process directly with no observation noise,\n", "and their task is to simulate the diffusion bridges between the observed values.\n", "Our examples uses observation noise with standard deviation σ_y > 0, as this is\n", "required for inference to work at this stage. =#\n", "const σ_y = 0.01\n", "\n", "#= The observation times ts and the corresponding observed values.\n", " It is assumed that the observation times are sorted increasingly. =#\n", "ts = [0.0, 1.0, 2.0, 3.0, 4.0, 5.0, 6.0, 7.0, 8.0, 9.0, 10.0]\n", "data_ys = [0.0, -0.009114765844660457, 0.029830170268001097,\n", "0.021755657377028295, 0.013934983652697013, 0.01910503829939146,\n", "0.049404191951179434, 0.03253467537333853, 0.07701445971172535,\n", "0.04618793002053633, 0.04001971135990065]\n", "N = length(ts)\n", "\n", "# The transition distribution of the OU process\n", "function trans(Yt, Δt)\n", " mean = Θ_1 / Θ_2 + (Yt - Θ_1 / Θ_2) * exp(- Θ_2 * Δt)\n", " variance = Θ_3 * Θ_3 / (2 * Θ_2) * (1 - exp(- 2 * Θ_2 * Δt))\n", " # Normal is initialized with mean and standard deviation, not variance!\n", " return Normal(mean, sqrt(variance))\n", "end\n", "\n", "# The generative model of the OU process\n", "@model ou begin\n", " # use TArray to support array copying when particles are duplicated!\n", " Y = TArray(Float64, N)\n", " Y[1] = 0.0\n", " for n = 2:N\n", " Δt = ts[n] - ts[n-1]\n", " @assume(Y[n] ~ trans(Y[n-1], Δt))\n", " @observe(data_ys[n] ~ Normal(Y[n], σ_y))\n", " end\n", " @predict Y\n", "end\n", "\n", "\n", "#= This concludes the example. The remainder of this script verifies that\n", "inference works properly by checking that (potentially weighted) samples from\n", "the posterior over Y[N] (the last element of the sequence Y) approximate the\n", "exact posterior over Y[N].\n", "\n", "The exact posterior over Y[N] is obtained by recalling that our OU process\n", "is a Gaussian process [Rasmussen and Williams, 2005] with constant mean function\n", "μ(t) = Θ_1 / Θ_2 and stationary covariance function\n", "\n", " k(r) = Θ_3^2 / (2 * Θ_2) * exp(- Θ_2 * |r|)\n", "\n", "which is also known as the exponential, Laplace, or Laplacian kernel. (See\n", "Appendix B.2.1 in [Rasmussen and Williams, 2005] for a derivation.) Hence exact\n", "posteriors can be obtained by performing standard GP regression, which involves\n", "inverting a kernel matrix A := (K + σ_y^2 I). =#\n", "\n", "k(x, y) = Θ_3^2 / (2 * Θ_2) * exp(- Θ_2 * abs(x - y))\n", "K = [k(x1, x2)::Float64 for x1=ts, x2=ts]\n", "noise_matrix = σ_y^2 * eye(N)\n", "noise_matrix[1, 1] = 0.0 # observed 0.0 at time t=0 with no noise\n", "A = K + noise_matrix\n", "Ainv = inv(A)\n", "\n", "# Compute true posterior over Y[N], the value of Y(t) at time t=ts[N]\n", "x_star = ts[N]\n", "k_star = [k(x, x_star)::Float64 for x=ts]\n", "f_star_mean = transpose(k_star) * (A \\ (data_ys - Θ_1 / Θ_2)) + Θ_1 / Θ_2\n", "f_star_variance = k(x_star, x_star) - transpose(k_star) * (A \\ k_star)\n", "f_star_std = sqrt(f_star_variance)\n", "f_star_posterior = Normal(f_star_mean[1], f_star_std[1])\n", "println(\"True posterior of last state: \", f_star_posterior)\n", "\n", "# KL-divergence between a Normal MLE-fitted to samples and true posterior\n", "function ou_divergence(samples, weights)\n", " fitted = fit_mle(Normal, samples, weights)\n", " return Turing.kl(fitted, f_star_posterior)\n", "end\n", "function ou_divergence(samples :: Vector{Float64})\n", " weights = fill(Float64(1), length(samples))\n", " return ou_divergence(samples, weights)\n", "end\n", "\n", "#= Compute the KL-divergence between moment-matched Normal from samples and\n", "exact posterior, as the number of samples increases =#\n", "num_samples = [10, 30, 100, 300, 1000, 3000, 10000]\n", "samples = sample(ou, SMC(maximum(num_samples)))\n", "samples_lasty = [res[N]::Float64 for res=samples[:Y]]\n", "kl_divergences = map(n -> ou_divergence(samples_lasty[1:n]), num_samples)\n", "\n", "# Plot the results using Gadfly\n", "using Gadfly\n", "plot(x=num_samples, y=kl_divergences,\n", " Scale.x_log10, Scale.y_log10,\n", " Guide.title(\"KL-divergence between exact and moment matched posterior\"),\n", " Guide.xlabel(\"number of samples\"), Guide.ylabel(\"KL-divergence\"))" ] }, { "cell_type": "code", "execution_count": null, "metadata": { "collapsed": true }, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Julia 0.4.6", "language": "julia", "name": "julia-0.4" }, "language_info": { "file_extension": ".jl", "mimetype": "application/julia", "name": "julia", "version": "0.4.6" } }, "nbformat": 4, "nbformat_minor": 0 }