{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# 7-MonteCarlo Sampling" ] }, { "cell_type": "markdown", "metadata": { "colab_type": "text", "id": "view-in-github" }, "source": [ "[![Open in Colab](https://colab.research.google.com/assets/colab-badge.svg)](https://colab.research.google.com/github/OpenJij/OpenJijTutorial/blob/master/source/ja/007-MonteCarloSampling.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "OpenJijはSimulated annealing (SA)を実行します。温度を一定に保てば、その温度でのカノニカル分布\n", "\n", "$$\n", "p(\\{\\sigma\\}) = \\frac{\\exp(-\\beta E(\\{\\sigma\\}))}{Z}, \\ Z = \\sum_{\\{\\sigma\\}}\\exp(-\\beta E(\\{\\sigma\\}))\n", "$$\n", "\n", "からスピン配列をサンプルすることが可能です。\n", "\n", "以下では全結合の強磁性イジングモデル\n", "\n", "$$\n", "E(\\{\\sigma\\}) = \\frac{J}{N} \\sum_{i" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# 結果の可視化\n", "plt.errorbar(temp_list, mag, yerr=mag_std)\n", "plt.plot(temp_list, mag)\n", "plt.xlabel('temperature', fontsize=15)\n", "plt.ylabel(r'$|m|$', fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "ある温度を境に値が大きく変わる、このような現象を**相転移**と呼びます。今回の模型では(システムサイズを無限大に近づけると)温度が1.0の時に相転移が起こることが、理論的に証明されています。 \n", "ただし、実際のモデルではどの程度の温度で相転移が起きるかを理論的に計算できない場合が多々あります。そのため、モンテカルロシミュレーションを用いて相転移に関する性質を数値的に研究する手法が多く取られています。" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Binder cumulant\n", "\n", "それでは、仮に相転移する温度を事前に知らなかったとして、数値計算からできるだけ正確に相転移する温度を求めてみましょう。上図を見ると、確かに磁化は温度が上がるにつれて0に近づいています。しかしどの温度が相転移点なのかまでははっきりしていません。これは、相転移現象は理論的にはサイズが無限大のシステムで起こる現象ですが、シミュレーションでは有限のサイズしか扱えないことにより理論との誤差が出てくることに起因します。これを**有限サイズ効果**と呼びます。 \n", "サイズが無限大のシステムの解析を数値的に行うことは、一見無理なことのように思われます。しかし、統計力学の数値計算分野に置いて、有限のシステムサイズから無限大のシステムサイズの情報を得る手法が開発されてきました。 \n", "そのうちの一つが**Binder cumulant**という量を用いる方法です。\n", "\n", "$$U_4 \\equiv \\frac{1}{2}\\left( 3- \\frac{\\langle m^4\\rangle}{\\langle m^2\\rangle^2} \\right)$$" ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [], "source": [ "# U_4の計算\n", "def u_4(states):\n", " m = np.array([np.mean(state) for state in states])\n", " return 0.5 * (3-np.mean(m**4)/(np.mean(m**2)**2))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "詳細は統計力学の教科書に譲りますが、この量は磁化が1に近づく強磁性では1に近づき、磁化が0に近づく常磁性では0になることが示せます。さらに相転移点ではシステムサイズに依存しない値を取ることが知られています。そのため、いくつかのシステムサイズで上述のような数値実験を行い、$U_4$のグラフが1点で交わる箇所が相転移点であると言えます。以下で実際の計算を行いましょう。" ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "# 系のサイズをリストで定義\n", "n_list = [40, 80, 120, 160]\n", "# 温度をリストで定義\n", "temp_list = np.linspace(0.5, 1.5, 30)\n", "\n", "# サンプラーを設定\n", "sampler = oj.SASampler(num_reads=300)\n", "\n", "u4_list_n = []\n", "for n in n_list:\n", " # インスタンス作成\n", " h, J = fully_connected(n)\n", " u4_temp = []\n", " for temp in temp_list:\n", " beta = 1.0/temp\n", " schedule = [[beta, 100 if temp < 0.9 else 300]]\n", " response = sampler.sample_ising(h, J, \n", " schedule=schedule, reinitialize_state=False,\n", " num_reads=100 if temp < 0.9 else 1000\n", " )\n", " u4_temp.append(u_4(response.states))\n", " u4_list_n.append(u4_temp)" ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAELCAYAAADHksFtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjMuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8vihELAAAACXBIWXMAAAsTAAALEwEAmpwYAABaNUlEQVR4nO3dd3gUxRvA8e/cXXrvkAQIEFpAOtI7UqUKCCIoUqRJExQ76g8VFEWk25AiVZAOIr1IB2mhJ0Co6b3c5eb3x8UQIIQkJDlC5vM89yS3O7v7TsR7b3dn3xFSShRFURTlUTTmDkBRFEV5uqlEoSiKomRJJQpFURQlSypRKIqiKFlSiUJRFEXJks7cAeQHd3d36efnZ+4wFEVRCpWjR4+GSSk9Hlz+TCYKPz8/jhw5Yu4wFEVRChUhxNXMlqtLT4qiKEqWVKJQFEVRsqQShaIoipKlZ/IehaIoyn/0ej0hISEkJSWZO5SnhrW1Nb6+vlhYWGSrvUoUiqI800JCQnBwcMDPzw8hhLnDMTspJeHh4YSEhFC6dOlsbWPWS09CiF+EEHeFEKcfsV4IIaYLIS4JIU4KIWoWdIyKohRuSUlJuLm5qSSRRgiBm5tbjs6wzH2PYj7QNov17YByaa/BwOwCiElRlGeMShL3y+nfw6yXnqSUu4UQflk06QwskKZa6AeEEM5CiOJSylv5Ec/Rnz4m9txZNEKL0GgRQoMQWjQa03uN0CLS1mk0OjRWNuisbdHZ2KK1ssbC2jbtZYeljR0WVrZorK0QOh1otQiN5uGfOt1974WVlemlMXcOVxRFMXna71H4ANczvA9JW/ZQohBCDMZ01kHJkiVzdbBrm9ZS/kzyfcuy83EtAUPaKzFXR36YQQt6ncBgocGQ/lNDqmXaT1srcHXC0t0TGy9vnIqVxM3XHzefsli6e6Cxs8ujSBRFMbdr164REBDAxIkTGTduHABHjx7l9ddfJzExkfbt2/P999/n25nT054oMut1pjMtSSnnAfMAateunavZmEJf+4KdV85hSE3BmJpMqiEZY6qeVGMyOmlAqzGgw4BGpKKVKeikHp1MwcKox0qTijUGLIQRC2lAYzSiMaYijEY0qUYwGhHpL5n2XprWGyVCpr1PFWhSQRhMPzWpAk2qBo1RizZVoE3VoNODzZ0kbK5EYB9/JT2ZxaS9AFIsNSQ72yKdHbCyssVaZ4O1zsZ0piIEiP9OP9P+xEIgbKzRubqhc3dDm/7TFZ27OzpXVzROTuoUXlHMYMyYMbRr1+6+ZUOHDmXevHnUq1eP9u3bs3nz5ofa5JWnPVGEACUyvPcFbubXwQZ3ag+0f2i5lJJkg5EkfSqJ+lQSUlJJSE4lLC6Z2zFJ3El7XY5JTv89LC7l3g4EoE17ZUGDEVdiKa6NpqpzElUckyhrE4ePRSxuMhLrpDCIuw1xdyElDoB4KbidasEdoyORwpG4FCuSkgTGOD3a2GSs4uLRJkmEBIHAWmuFjdYaa6011lorrDVWWOmsEAjk3bskHj9BamQkGI0PB6jToXN1RevuhoWHJzovL3Senui8PLHw8kp/r3V2VglFUTIIDg6mXbt2NGrUiP379+Pj48OaNWuwsbF57LZ//vknZcqUwS7DVYJbt24RExND/fr1AejXrx9//vlnkU0Ua4ERQoilQF0gOr/uT2RFCIG1hRZrCy3O2dxGn2okNNaUSOKSDGg1AiFAK0Ta76afGgGatN9TjZLLoXGcux1L4K0Ytt2K5XbMvZEJ7vZWVCruQCV/R54vpqGxRwJ2sdcpG3WVspFXIeoqRF6FqGvgZroIliAEwXYuXCleiStOXlzWCa4khnI99jqpMtXUPwQ+9j5423tja1EKe60tLkk6nBMEjvFGHOKM2MbpsYlJxiomCYuoeCxv34TTp0kND3/472VlZUognp7YVKmC6xtvYOHl+aT/GRTliX267gxnb8Y8vmEOBHg78knHyo9td/HiRZYsWcKPP/5Iz549+eOPP7h16xaLFy9+qG2TJk2YPn068fHxTJ48ma1bt/LNN9+kr79x4wa+vr7p7319fblx40bedCgTZk0UQoglQDPAXQgRAnwCWABIKecAGzF9xb8EJAD9zRNpzlloNXg72+Dt/PhvDBlV8XGic4b3EfEpnLsVQ2Ba8jh3O4b5+4OZZzBib6WjdeWSdKpWj4Z13bHQpl2EktJ01hF1FdvwSwRc3UdA0B64uM+03s6TFL+GXPOuwmUnL64YYrkcfYW7CXe5GXeTeH08CfoE4vXxpBhTwAbT64GakpYaS0ra+FPB6EFpvTO+iTZ4xmtxijViFZmI8c5dIhYvJnLJEpx7vYzbwIFYeKqEoRRNpUuXpnr16gDUqlWL4OBgPvzwQ8aPH//IbT755BPGjBmDvb39fctN43vul59n8eYe9dT7MeslMLyAwnkqudpZ0sDfnQb+7unL9KlGDl6JYO2/N9h0+jarjt3A1c6SDs8Vp1N1b2qVdEHj4AUOXlDieaj+imnDyGAI2gNBu7EM2o3/mdX4Azj6Qukm4NcafGqDe3lIG3WlN+rTk8Z/r9iUWG7F3+J67HWux17nYux1tutPk2iRCM6AM4gSAi87L1p2b0XvAxZELv6dqGXLcenVC7dBA9G5u6MoBS073/zzi5WVVfrvWq2WxMREvv766yzPKA4ePMjKlSt55513iIqKQqPRYG1tzUsvvURISEh6+5CQELy9vfMt9qf90pOSCQuthkbl3GlUzp3Pu1Rh5/lQ1v57kxVHr7PwwFW8nazpWM2bTtW9CSjueO+bhouf6VWzr+msI+wiBO+GoN1wYTP8+7upnZUjeNcA39pY+NTGyacWTg5Z/yOUUhKeFE5IbEh6AgmOCWZZ8Fa21/bgqx6TKbZiLxELFxK5bBkur7yC24A30Lm55evfSlGeZuPHj8/yjGLPnj3pv0+cOBF7e3tGjBgBgIODAwcOHKBu3bosWLCAt956K9/iVImikLPSaWlTuRhtKhcjLtnA32fvsPbfm/y8N4i5u6/g72nPa/VL0aN2CawtMtxNFwI8yptedQaabl6HX4SQI3DjiOnn3mmQdh8DpxLgUwt8a4NvHdPvWosMuxO427jjbuNOdc/q6ctPB5zmnd3v0P/M+wzuNZj+g9cSNXceEfPnE7lkCa59XsF1wAB0Li4F8wdTlGfE7Nmz04fHtmvXLt9uZAOIzK51FXa1a9eWRX3iosj4FDaevsWKIyGcuB6Fu70l/RuW5tV6pXCyyV4hMFIS4PbJDMnjKERfM62zdobybaBiByjbEqzsH7mbeH08Xxz8grWX11LdozpfNfkK97sphM2aRcyGDQgbG1xffRX3Nwer5z+UPBcYGEilSpXMHcZTJ7O/ixDiqJSy9oNtVaJ4xkkpORQUwexdl9l5PhR7Kx196pVkQMPSeDpa53yHcXfh2j9wfjNc2ASJkaC1grLNoUJ7qNAO7DO/Yb3xykY+P/A5AB/X/5h2pduRfOmSKWFs2oxlmTL4fj8NK3//J+myotxHJYrMqUShEkWmztyMZs6uK2w4eROdVkP3Wr4MblwGP/dcfotPNcD1A3BuA5xbbxqWi4ASdU1nGhU7gFvZ+za5EXeDd3e/y7+h/9KpbCfer/s+dhZ2xB84wI23x2FMSKD4Z5/i1LHjk3dYUVCJ4lFUolCJIkvBYfHM23OFlUdCMBiNtH+uOEOalqWKj1Pudyol3DkN5zaaksbtk6blpRpCt3ngdG/Mt8FoYO7Jucw7OQ9fe18mN5lMFfcq6O/c5cbbY0k8chTnXi/j9d57aDKMFFGU3FCJInMqUahEkS13Y5L4eV8Qiw9cIy7ZQM2SztQp7Uqtki7ULOWCu/0TfEhHXYOza2DnV6Czgpd+Nl2eyuDonaNM2DOBsIQwRtQYQf8q/RGpRkKnTSP8p5+xDgjAZ/r3WGZ4sEhRckolisypRKESRY5EJ+pZdOAqW8/e4czNaPSppn8TpdxsqZmWNGqVdKFCMQe0mhw+1BN2EZa9CqHnocUH0Ojt9Gc0AKKTo/n0n0/ZenUr/Sv3Z2ztsQDEbt/OzQnvAeD91Zc4tGiRN51VihyVKDKnEoVKFLmWpE/l9I1ojl2L5OjVSI5ejSIszlRR185SS7USztQu5UK3mr7Zv7eRHAfrRsHplVCuDXSbCzb3hsNKKZl0cBLLzi/jfw3/R2d/07PpKSEh3Bg5iqSzZ3EbOACP0aNNJdsVJQdUosicShQqUeQZKSUhkYnpiePYtUgCb8UipaTdc8UZmt17G1LC4Z9g83vg6A09F4B39fTVeqOeoVuHcuzuMX5p80v6sxjG5GTufPklUUuXYVO7Fj5Tv1V1o5QcKcyJQq/XM3DgQI4dO4bBYKBfv368957pTPtJy4znJFEgpXzmXrVq1ZJK/rkTnSi/3Bgoq3y8WZZ6d73s8+MBuedCqDQajY/f+NohKadWkvIzDymP/nbfqqikKNn+j/ayydIm8kbsjfvXrV0rA6vXkOcbNJRxBw7mZXeUZ9zZs2fNHUKuLV68WL788stSSinj4+NlqVKlZFBQkJRSyjp16sj9+/dLo9Eo27ZtKzdu3JijfWf2dwGOyEw+U9U0akqOeTpaM6FdRfa914IJ7Spy/k4sr/58kI4z9rL+5E1SjVmcpZaoA2/uhpL1YO1bsGY46E2Vbp2snPih5Q/oU/WM3D6SBH1C+mZOHTtSesVytM7OXH/zTRJPn8nvbipKngkODqZSpUoMGjSIypUr07p1axITHz/NmRCC+Ph4DAYDiYmJWFpa4ujoeF+ZcSFEepnx/KIu+Cq55mhtwZCmZenf0I/Vx24wb/cVRvx+nFJu5xnUuAzda/neXzbkP3bu0Hc17PgC9nwDt/6FngvBtTRlnMowpekUhm8bzvt73+fbZt+iEabvM1b+/pT6bT5BPXsSMmwYfiuWY+HlVcC9Vgq1TRPg9qm83Wex56DdV49tlpsy4927d2fNmjUUL16chIQEvvvuO1xdXTly5EjRKTOuPBusdFp6PV+SHrVLsPXsbWbvusKHf55m2t8XeKdNRXrWKfHwRhottPzIVDdq9WCY28T0zIVLKRo5l2JcifZMubaeGUe+ZWSdcemb6dzdKTF7Dld79yZk2HBKLVqIJhuTvyiKueWmzPihQ4fQarXcvHmTyMhIGjduTKtWrYpWmXHl2aLVCNpWKU6bysU4cCWCb/46z3urT1GzlDP+ng6Zb1ShLQzeBds+Mw2hDd4DKXG8Clx2d+XHs79R9p+5dLD2BZdS4FwK60qd8J76DSHDhnPz3Qn4TPvONMWrojxONr7555fclBn//fffadu2LRYWFnh6etKwYUOOHDlC48aNC7TMuPq/S8lzQgjql3VjXt9a2FpombQhMOsNXEtDj19h2H54LwTGX0EM2s4HTSZTy7oYHzvZcsra2nSJ6p+ZsLALDvVr4vnOO8T+9RehP/xQMB1TlDw2fvx4Tpw48dBr+vTpAJQsWZLt27cjpSQ+Pp4DBw5QsWJFihcvnl5mXErJggUL6Ny582OOlnsqUSj5xs3eihEt/NlxPpTdF0Kzt5EQYOcGPrWwqNqT7zovx8O+OKOsk7kzYDMM/BuSY+Dwz7i+/hrOPboTPnsO0evW5W9nFMUMhg8fTlxcHFWqVKFOnTr079+fqlWrAqYy4wMHDsTf35+yZcuqMuM5pZ6jeHokG1J54dvdWFto2DiyMTptzr+bXIy8yKsbX8XPyY/5bedjs6SPqZbU6FNIqeXawEEknjhByd/mY1ujRj70QinMCvNzFPkpJ89RqDMKJV9Z6bS8374iF+7EseTw9Vzto5xLOaY0mUJgeCAf7fsI2WgMxIfC8UUIS0t8vp+GrngxQka8hT4fR34oSlGlEoWS79pULkbd0q58t/UC0Yn6XO2jaYmmjKk1hi3BW5gXc9ZUynzfdEjVo3NxocTs2ciUFK4PGUpqXHwe90BRijaVKJR8J4TgoxcDiExIYcb2i7nez+uVX6edXzvmnJzD9dr9TLPtnVoJgFXapEfJV65w8+23kampeRW+ohR5KlEoBaKKjxPda/oyf38wwWG5+8YvhGBcnXFYaCyYHvUveFaGvd+Z5vsG7Bo0oNiHHxC3axd3v/4mL8NXlCJNJQqlwIxvUwELrYYvNz1muGwWPG096RfQj83BmzlV82UIOw/nN6avd+ndG5dXXyVi/nwiV6zIi7AVpchTiUIpMJ6O1gxrVpYtZ+6w/3JYrvfTv0p/XK1dmRp5DOlSCvZMNVWnTeM14V3sGjXi9qefkXjyZF6ErihFmkoUSoEa2LgMPs42/G99YNbFA7NgZ2HHsGrDOHr3GLuqdoSbxyBoV/p6odPh89236Fxduf2/Sci0S1OKUtiEh4fTvHlz7O3tGTFiRPryhIQEOnToQMWKFalcuTITJkxIX5ecnMzLL7+Mv78/devWJTg4+InjUIlCKVDWFlrebVeRs7diWHk0d8NlAbqV74afox/fRZ3CYF8M9nx733qtgwOeb48l6eRJoteufdKwFcUsrK2t+fzzz/nmm4fvuY0bN45z585x/Phx9u3bx6ZNmwD4+eefcXFx4dKlS4wZM4Z33333ieNQiUIpcB2rFqdmSWe+3nKBuGRDrvZhobFgdM3RXIkJ4s+AlqYzipCj97Vx7NgR66pVCZ36LcZ4NWRWMZ/clhm3s7OjUaNGWFtb37fc1taW5s1Nc9BbWlpSs2bN9NpPa9as4bXXXgOge/fubNu2LdMigjmhigIqBe6/4bJdZ+1n9s5LjG9TMVf7aVGyBTU8azAz5iztbZyx3fst9LpXYE1oNBR7/z2Ce/UmbN6PeI4ZnUc9UAqryYcmcy7iXJ7us6JrRd59/vHf2nNTZjw7oqKiWLduHaNGjQLgxo0blChhqtis0+lwcnIiPDwcd3f3HPTqfuqMQjGLGiVd6FLdmx/3BHE9IuHxG2RCCMHYWmMJSwrntwoN4dx6uHv/h4BN9eo4dupIxK+/kpKh2qaiFLTMyow/rijg4xgMBnr37s3IkSMpU6YMQL6UIFdnFIrZvNO2IpvP3Gby5nPMeKVmrvZR3bM6L5R6gV9D9tLDyg73vd9Bt7n3tfF8+21it/7N3Slf4zv9+7wIXSmksvPNP7/kpsz44wwePJhy5coxevTo9GW+vr5cv34dX19fDAYD0dHRuLq6PlHsZj+jEEK0FUKcF0JcEkJMyGS9kxBinRDiXyHEGSFEf3PEqeQ9b2cbBjcpy/qTtzgSHJHr/YyqOQq9Uc+csjXh1AqIvHrfegsvL9wHDyL2r7+IP3joScNWlDzzJGcUH374IdHR0UybNu2+5Z06deK3334DYOXKlbRo0eKJzyjMmiiEEFpgJtAOCAB6CyECHmg2HDgrpawGNAOmCiEsCzRQJd8MaVoGL0crPl9/FmMuh8uWcixFjwo9WJl0nSALC9j/8PwUrv37Y+HtzZ0vv1TlPZRCxc/Pj7FjxzJ//nx8fX05e/YsISEhTJo0ibNnz1KzZk2qV6/OTz/9BMCAAQMIDw/H39+fb7/9lq++evLJmsxaZlwIUR+YKKVsk/b+PQAp5ZcZ2rwHlMCUMPyArUB5KeUjB8erMuOFyx9HQ3h7xb/UKuXCkKZlaVnRE40mZ9+AwhPD6bC6A3Wx5vsrgTD6FNh73tcmZvNmboweQ7FPP8Xl5Z552QXlKabKjGeuMJUZ9wEyDqYPSVuW0QygEnATOAWMyixJCCEGCyGOCCGOhIZmc5Ic5anQraYPn3euzO3oJAYtOMIL3+1i2eFrJBuy/83fzcaNN6q8wXZ9OMd0wIFZD7VxaNMGm9q1CJ02jdSYmDzsgaI828ydKDL72vjgKU4b4ATgDVQHZgghHB/aSMp5UsraUsraHh4eeR2nko+EEPSt78eu8c34vld1rHRa3v3jFI0n72D2zsvZLk3eN6AvnjaeTPUpjTz8MyRGPXScYu+/T2pUFGGzZudDTxTl2WTuRBGC6bLSf3wxnTlk1B9YJU0uAUFA7gbeK081nVZD5+o+bBjZiIUDnqe8lwOTN5+j4VfbmbThLLeis35AyUZnw4gaIzhpjOdvrR4O//RQG+uAAJy7v0TEokUkXwnKr64oyjPF3IniMFBOCFE67QZ1L+DBegvXgJYAQggvoAJwpUCjVAqUEILG5TxYNLAu699qRIuKnvy8N4jGk3fw9vJ/uRIa98htO5XthL+zP9O8vNEfmA0pDz+j4TFqFBpra+5Onpyf3VCUZ4ZZE4WU0gCMALYAgcByKeUZIcQQIcSQtGafAw2EEKeAbcC7Usrclx5VCpUqPk5M712DXeOb82q9Umw4dZOus/ZzLTzzh/S0Gi1jao3hGnpWaJPgxMNj1HXu7rgPHUrcrl3E7dmT311QlELPrKOe8osa9fTsCg6Lp9OMvfi42LJqaANsLLUPtZFSMvCvAVy8dYQN0eDw1nHQ3N9OpqRwpWMn0Gops+ZPhIVFQXVBKWBq1FPmCtOoJ0XJET93O77vXYNzt2OYsOrkI8sVvF17HDFCMM4yAf2ZVQ+3sbTE8913SblyhcglSwoidEXJsUeVGQdISUlh8ODBlC9fnooVK/LHH38Aqsy4ogDQvIInb79QnjUnbvLLvuBM2wS4BTCxwSfst7Xh/cNfkZr6cJVa++bNsGvYkNAZMzFERuZv0IqSC1mVGZ80aRKenp5cuHCBs2fP0rRpU0CVGVeUdMOa+dM6wIsvNgY+cra8LuW68bZnYzZrU/hy+6iHzj6EEHi9NwFjfDxhPzz8NLei5JW8LjMO8Msvv/Dee+8BoNFo0qvDqjLjipJGoxFM7VmNLjP38dbvx1n3ViO8nW0eavd6q6lE/FiDX2/uxvXf2QyrPuy+9Vb+/jh3707UipW4Dx+Ozs2toLqgmMHtL74gOTBvy4xbVapIsffff2y7vCwzHhUVBcBHH33Ezp07KVu2LDNmzMDLy0uVGVeUjBysLZjbtzbJBiNDFh0lSZ/Jk9wWNoyp9DpdY+OY/e9sfg/8/aEmrq+9htTriVy6tACiVoqqvCwzbjAYCAkJoWHDhhw7doz69eszbtw4QJUZV5SH+Hva823PagxeeJSP/jzNlO5VH/qfQjw/iI/3TSPKpRRfHfoKZytn2pdpn77eqkxp7Bo3JnLpUtwHDUJYqpqTz6rsfPPPL3lZZtzNzQ1bW1u6du0KQI8ePfj555+BZ7TMuKI8qdaVizGyhT8rjoaw6OC1hxvYuaGr8SpTrpymplsVPtj7Aftu7LuviWu/vqSGhhGzeXMBRa0ouS8zLoSgY8eO7Ny5E4Bt27YREGAqvP3MlRlXlLwyulV5mlfw4LN1ZzKf26L+cKyNqfxg4Ye/iz9jdo7h39B/01fbNWyIZenSRCxY+MQ3/hQlL2VWZhxg8uTJTJw4kapVq7Jw4UKmTp0KPINlxvOLeuCuaIpO0NNp5l4SUlJZ/1YjvBwfGCmy/DW4vIOwYXvot20oMSkx/Nb2N8o6lwUg4vffufPZ55T6/Xdsa9YwQw+U/KAeuMuceuBOKZKcbC2Y17c28ckGhi0+RorhgWr0DUdCcjTuZ9Yx94W5WGgsGLx1MDfjTHUonTt3RuPgQMTCBWaIXlGeXipRKM+UCsUcmNK9KkevRvLZ+jP3r/SpBaUawYHZlLAtxpxWc0jUJ/Lm1jeJSIpAY2eHc/fuxP61Ff2tW+bpgKI8hVSiUJ45L1b15s2mZVh04BqbTz/wgd9wJMSEwOlVVHCtwIyWM7gVf4t3dr8DgEufPiAlkb+rsh7PkmfxEvuTyOnfQyUK5Zk0vnUFKhZzYNLGwPufr/B/ATwqwv7pICU1vWoyvPpwDt46yMXIi1j6+uDQsgVRy5djzMaTs8rTz9ramvDwcJUs0kgpCQ8Pz/Rp70dRz1EozySdVsMHHSrR9+dDzN8fzJCmphvWaDTQ4C1YMxwubwf/lnTx78KM4zNYcWEF79d9H5e+fYnd+jfR69bh0lPNrV3Y+fr6EhISgpoi+R5ra2t8fX2z3V6NelKeaQPmH+ZgUAQ7xzfD3T7tgSdDMkyrCp4Vod8aACbsmcCu67vY1mMbNjobgrp2g1QDpdeufeIx6IpSWKhRT0qR9H6HSiTpU/l264V7C3VWUG8IXNkJt0zPUvQs35M4fRxbgrcghMC1b1+SL14i4cAB8wSuKE8RlSiUZ1pZD3terVeKpYeuce52zL0VtfqDpT3sN1WNreFZg7JOZVl+fjkAji92QOvqSsSCheYIW1GeKipRKM+80a3K4WBtwf/WB967oWnjDDVfg9OrIOoaQgh6VOjB6fDTnA0/i8bKCueXexK3cycp1zIpC6IoRYhKFMozz9nWktGtyrH3Uhjbz929t6LeUNPPA7MB6Fi2I9Zaa1ZcWAGAS6/eoNUSsWhRQYesKE8VlSiUIuHVeqUo42HHpI2B6FPTnth2LgFVXoKjv0FiJI6WjrQt3ZaNVzYSr4/HwssTx7Ztif5jFalxcebtgKKYkUoUSpFgodXwQftKXAmNZ9GBq/dWNBwJ+ng48gsAPcr3IMGQwIYrGwBTVVljfDzRq1abI2xFeSqoRKEUGS0qetLI351pf18kKiHFtLDYc1C2pemmdnw4z7k/R0XXiiw/vxwpJTZVq2JTrRoRixchjcasD6AozyiVKJQiQwjBhy9WIjZJz7S/L95b0fp/kBQD2yaabmqX78H5yPOcCjsFgEu/vuivXiNu1y4zRa4o5qUShVKkVCzmyMt1SrLowFUuh6bdd/AKgPrD4NgCuH6IDmU6YKuzTb+p7di6NTovLyIXqqGyStGkEoVS5LzdujzWFlq+2BB4b2HTCeDgDRvGYqexon2Z9mwO2kxMSgzCwgKX3r2J3/8PyRcvPnrHivKMUolCKXLc7a0Y0cKfbefusvdimGmhlT20/RJun4LDP9GzfE+SUpNYd3kdAM4v90RYWRGxUA2VVYoelSiUIql/Qz9KuNrwvw1nSTWmPYQX0Nl0Y3v7/6hk6UIVtyqsOL8CKSU6FxccO75I9Nq1GCIjzRu8ohQwlSiUIslKp+W9dpU4dzuWZYevmxYKAe2/htQU2PIBPSv05HL0ZY7fPQ6Aa79+yJQU7n79jRkjV5SCpxKFUmS1q1KM5/1cmfrXeWKS9KaFbmWh0Rg4vZI22GFvYc/yC6b6T9bly+M2eBDRq1YR89dfZoxcUQqWShRKkfXfcNnw+BS+3BhIsiFtgqNGo8HFD9vN79OxdAe2Bm8lMsl0uclj2DCsAwK4/fEnGNT8BkoRYfZEIYRoK4Q4L4S4JISY8Ig2zYQQJ4QQZ4QQajC7kmeq+jrzegM/lhy6TotvdrH88HUMGito/w2EX6RHQhIpxhTWXl4LgLC0xPvrKRgTE7n5wQdq1jSlSDBrohBCaIGZQDsgAOgthAh4oI0zMAvoJKWsDPQo6DiVZ9vETpVZPLAu7g5WvPPHSVpP282GxCrIih0p989P1HANYMWFFelJwapsWTzHjyd+9x6ili41c/SKkv/MfUbxPHBJSnlFSpkCLAU6P9DmFWCVlPIagJTyLoqSxxr6u/PnsAbM7VsLrRAM//0Y/W93IxVB9+horsZc5dDtQ+ntXV7pjV3DhtyZPIXkoCAzRq4o+c/cicIHuJ7hfUjasozKAy5CiJ1CiKNCiH6Z7UgIMVgIcUQIcUTNjavkhhCCNpWLsXl0E77tWY3LKc5MSexM68sHsNfYpE9qBCA0Gop/8QXCyoqb77yL1OvNGLmi5C9zJ4rMJiN+8KKvDqgFdADaAB8JIco/tJGU86SUtaWUtT08PPI+UqXI0GoE3Wr6sm1sM0q2f5ub0ocXImLYGvw3+zKcPVh4eVL804kknTpF2Ow5ZoxYUfKXuRNFCFAiw3tf4GYmbTZLKeOllGHAbqBaAcWnFGGWOg19GpbDp89s+seFIoWR/itnMn/fvWTh2LYtTp07ETZ3LoknTpgvWEXJR+ZOFIeBckKI0kIIS6AXsPaBNmuAxkIInRDCFqgLBKIoBcSqXBNKV+pOnaRkXDwOMnHdaX7NkCy8PvwQnZcnN959F2N8vBkjVZT8YdZEIaU0ACOALZg+/JdLKc8IIYYIIYaktQkENgMngUPAT1LK0+aKWSmiWn9Oz4RUEkQUz1e6y6frzvLLXlOy0Do44P3VV+ivXefOlK/NHKii5D3xLI4Dr127tjxy5Ii5w1CeMfqDc2h7+nvK2BbDwjCZjWdC+bBDJQY2LgPAna+/JuLnX/CdPQuH5s3NHK2i5JwQ4qiUsvaDy8196UlRCg2LOoPo7ViJAylhjOIzugXY878Ngfy05woAHqNGYVWhArc+/AhDeLiZo1WUvKMShaJkl0ZLj06/YiN0LIoJZGr027xeIZX/bQjkx91X0Fha4j1lCsaYGG59/Il6alt5ZqhEoSg54GTlROfy3dng4Eh4UgSf3BnB22VvMGljIPN2X8a6Qnk8xowhbts2ov/4w9zhKkqeUIlCUXLo1YBXMchUljZ6A+Hoy4ibE/i6xH6+2BjInF2XcX39NWxq1+Lud9OQKSnmDldRnphKFIqSQ6UcS9GsRDOWXd1C0mtrERXa0SN0Br97LebbTaeYszsIt4EDSQ0PJ3b7dnOHqyhPTCUKRcmFfgH9iEqOYt2NHdBzITR5hwbRG9nk8jU/bT7IQmNxdN7FiVy2zNyhKsoTU4lCUXKhllctAtwCWHh2IUYBtPgAuv9KGf1lttpPZP3Wv7lUuwUJ/xwg5epVc4erKE9EJQpFyQUhBP0C+hEUHcTeG3tNC6t0Q7yxGRdbHautP+MvQxxSoyVy+fKsd6YoT7k8TRRCCLe83J+iPM1a+7XG09aTBWcX3FvoXR0xeCc6r4p85vgbp3wrEbVqNUZ1U1spxPL6jOJwHu9PUZ5aFhoL+lTqw8FbBzkfcf7eCntPtJ2m4UACCSXBGBlJ7Nat5gtUUZ5QrhOFEMIohEhN+2kUQqQCfhl+V5Rn3kvlXsJGZ3P/WQWAdw0o25IepQ5zx9aFK78sMk+AipIHnuSMYg6wEPCQUmqklFrgaobfFeWZ52TlRFf/rmwM2khowgMTZjUZh31qJAmVvLA+c4Jbp86ZJ0hFeUK5ThRSymHAT8CfQojX/lucJ1EpSiHyaqVXSTWmsuTckvtXlGoAJevTvMI5DELD1ilzVVkPpVB6onsUUsq9QEtMc0psBSzzJCpFKURKOJagRckWLL+wnERD4v0rG4/DxniT1EqlqHBiN7/vuWSeIBXlCTw2UQghZqXND9FQCOHw4HopZYqU8kNgLPC//AhSUZ52/QL6EZ0czbrL6+5f4d8SilfDv+wNHPUJ7P5pOVdC48wTpKLkUnbOKDoCszBNQRolhAgSQqwRQnwuhOj235BYKeUpKaWaOFgpkmp41qCKWxXTA3jSeG+FEND4bexsg9F6utA66B/GLDuBPtX46J0pylPmsYlCSlkCcAWaA6OArYAXMBpYCdwSQiwTQvjmY5yK8lQTQtCvcj+CY4LZE7Ln/pUVOyI8yuNaLomA0MtEBF5gxnZ1CUopPLJ1j0JKGSWl3C2lnCGlHCylrCeldAAqAROA+sBBIYR3fgarKE+zVqVaUcyu2MNDZTUaaDQWZ48g0Gl5K/40M3Zc4vi1SPMEqig59KQ3s89LKb8FqgMGYGIexKQohZKFxoI+Fftw6PYhAsMD71/5XHd0xXxxKGtJ1bP7KGGrYcyyE8QnG8wTrKLkQJ48mS2ljACmA+3yYn+KUlh1K98NW53tw2cVWgtoOAoX7+vImBimeoVzNSKB/20IzHxHivIUycsSHlcB9zzcn6IUOo6WjnQr143NQZu5E3/n/pXVX8W2jAuWLha47djI4MZlWHLoGtsC728Xl2zg4p1Ydp6/y+8Hr/HNlvOMXXaCl+f+w9S/zqtnMZQCp3tcAyHERuAk8G/az3NSysxKdNQFrudteIpS+LxS6RV+P/c7i88tZmytsfdWWFgjGozA+ehX3D12jBEfa9h90ZHxK09Ss6QzN6KSuBmVSHSi/r79aTWCYo7WOFjr+GH7Jbydbej9fMkC7pVSlD02UQAOwJuAE6Ynr1OEEIGYksZZIAKoCQwAvsynOBWl0CjhUII2fm1YcGYBNT1r0qxEs3sra7+BU8WphJ4SxK9cyfeDRjJyyXFCIhPxcbahdikXvJ1t8Ha2xsfZBm9nGzwdrNBpNRiNktd+PcQna89QzdeZAG9Hs/VRKVpEdk9jhRAlgaoZXs8B5QEtkAz8AoyRUpq9nnLt2rXlkSNHzB2GUoTF6+MZuGUgF6MuMrvVbOoUq3Nv5c6vuPHFbOIi3Cm3Zy8aG5ts7zcsLpkO0/dgZ6lj7VuNsLfKznc9RckeIcRRKWXtB5dn+x6FlPKalHK9lPILKWUvKWVlwAbwAeyklMOfhiShKE8DOws7ZrWahY+9D29tf4uz4WfvrXx+MM4VJca4BGI2bc7Rft3trZjeqwbB4fG8v+qUul+hFIgnHR5rkFLeklKqx0wV5QEu1i7MfWEujpaODP17KEHRQaYVtq7Yvvgalo4Gohb/luP91i3jxtgXyrP235ssPaxuCyr5T02Fqij5qJhdMea9MA+AwVsHczv+NgCi/gic/ZNJPHOepPPns9pFpoY186dxOXc+WXuGszdj8jRmRXmQShSKks/8nPyY02oOcSlxDN46mMikSHDwwrlrF4RGErXg5xzvU6MRfPdydZxtLBjx+zHi1IN7Sj5SiUJRCkAlt0rMaDmDm3E3Gfr3UOJS4tC+MA6HkslEb9iEMSEhx/t0t7diem91v0LJfypRKEoBqeVVi2+bfcv5iPOM2jGKZAcvXNrWx5hkIPK3H3O1z3rqfoVSAMyeKIQQbYUQ54UQl4QQE7JoVydtju7uBRmfouSlJr5N+LzR5xy6fYjxu8Zj+con2HnrCZszB8O13E2V+t/9iolrzxB4S92vUPKeWROFEEILzMRUIyoA6C2ECHhEu8nAloKNUFHy3otlXuS9599jx/UdfHJhAR4fT8KYIgl7uztEh+R4f//dr3CysWD4YnW/Qsl75j6jeB64JKW8kvYMxlKgcybt3gL+AO4WZHCKkl9eqfQKw6oPY+3ltUy3Dca5fQsiTxtI/q4tRAbneH8Z71d8sFrdr1DylrkThQ/314cKSVuWTgjhA3QFspw9TwgxWAhxRAhxJDQ0NM8DVZS8NqTqEPpU6sOiwEUsaeOOsLbh7r5k+LU9hOV8YqP/7lesOXGTZep+hZKHzJ0oRCbLHvwqNA149xGFCO9tJOU8KWVtKWVtDw+PvIpPUfKNEIJ36rxDv4B+/HprNQdfKEHcdR3x11Ngfnu4m/N7Fhmfr1h1LESdWSh5wtyJIgQokeG9L3DzgTa1gaVCiGCgOzBLCNGlQKJTlHymERrG1xnPO3XeYXr5IKKdLbl5qRxSClOyuH0qZ/tLu19RqbgjY5f/S/c5/3D6RnQ+Ra8UFeZOFIeBckKI0kIIS6AXsDZjAyllaSmln5TSD9Mc3cOklH8WeKSKko/6BvTly1ZTWdwUDBeDCXYbBDprmP8i3DiWo32521uxamgDpnSvSnBYPB1n7OWD1aeIjFel2JTcMWuikFIagBGYRjMFAsullGeEEEOEEEPMGZuiFLTWfq3pP/pngry13JrzK2fafQ3WTrCgM1w7mKN9aTSCnrVLsH1cM15v4MfSw9dpPnUniw5cJdWoLkcpOZPtMuOFiSozrhRml3atQ//mO/zR1JKmb0+k8ebPIfY29FkOfo1ytc/zt2P5ZO1pDlyJoLK3I592qkxtP9c8jlwp7J64zLiiKAXDv2lHLFs05cV/DHz090RWNX8LnHxhUXc4Oh9unYTkuBzts0IxB5YMqscPvWsQEZ9C9zn/MHbZCe7GJOVPJ5RnijqjUJSnUMrVq1zu8CKn67jzWbMwhlTqx7CjfyLunL7XyN4LXMuCaxlwK2P6+d97K/tH7jshxcDMHZf4cXcQljoN3Wr64ONsg5ejNZ6OVhRztMbL0Ro7NSlSkfOoMwqVKBTlKXXny6+IWLCAjZ+0ZH7yTrqU7cTHZbpjEXkVIi5DxBUIv2L6GXf7/o2dS0LPheBd/ZH7Dw6LZ9LGQP65HJ7p09z2Vjq8HK3wSkscHg5WWOk06DQadFqBhVag1Wiw0ArTMo1ApxXotBpqlHCmhKttHv9FlPymEoWiFDKpUVFcatMW64AANoysxeyTc2jn144pTac83Dg5DiKDIDwtgRycYzrjGLQDtI8/M4hLNnAnJinDK5nb0Uncjb33e2hcMimG7M1R5mJrwfqRjfFxzv40r4r5qUShKIVQxIIF3PniS0rMncN8x1PMPTmXhe0WUt2zetYbnlkNK16HNl9C/WF5Fo+UEqMEfaqRVKPEkCrRG02//7fsbmwyb/x6mNIedix/sz7WFto8O76Sv9TNbEUphFx69cKiVEnuTPma/hX74WrtyowTMx6/YUAX8H8BdkyC6Bt5Fo8QAq1GYG2hxc5Kh5OtBe72pstTvi62lHKzo46fK1N7VuNkSDSfrjuTZ8dWzEclCkV5iglLS7zGjyfl8mVS/tzAgCoDOHjrIIduHXrMhgI6fANGA2x+t2CCzaB15WIMb16WJYeus+zwtQI/vpK3VKJQlKecfcuW2NauTegPM+jh0x5PW09mnJjx+DpOLn7Q9B0IXAfnNxdIrBmNfaECjfzd+WjNGU6FqDIihZlKFIrylBNC4Pnuu6RGRBD703zerPomx+8eZ9/NfY/fuP5b4FERNo6DlPj8DzYDrUYwvXcNPOytGLLoqCohUoipRKEohYDNc1Vw6tqV8F/n0y6xHD72Pvxw/IfHn1XoLOHFaRB9HXZNLpBYM3K1s2RWn5qExiYzatkJVT6kkFKJQlEKCa/3JqDz8ODuhPcZWuENzoafZfu17Y/fsFR9qNEX/pkJdwr+5nK1Es582rkyuy+E8v3fFwr8+MqTU4lCUQoJraMj3l9+QUpwMHVWBeLn6MeMEzNINWY5VYvJC5+ZCgyuGw3G7D0LkZd61SlBz9q+TN9+iW2Bdwr8+MqTUYlCUQoRu/r1cX2tH9G/L+Xt1JZcirrEluBsTCVv6wqt/wchh+D4gvwP9AFCCD7rXIUqPo6MXnaC4LCCvV+iPBmVKBSlkPEYMwbLsmXx/n4VVS1LM+vfWRiMD5fgeEi13uDXGLZ+AnEFP12wtYWW2X1qodUIhiw6SmJKNs6ElKeCShSKUshorK3xnjKZ1MhI3t7pwNXoYNZdXvf4DYWADt+aRj/99UH+B5qJEq62fN+rBufvxPL+6lNqqtZCQiUKRSmEbCpXxmPECGx2HePlaz7M+XcOKanZGH7qUR4ajYaTy+DKzvwOM1NNy3swtlV5Vh+/wcIDV80Sg5IzKlEoSiHlNnAANjVq0HVNKMm3brDq4qrsbdj4bXApDevHgt4881EMb+5Py4qefLbuLDN3XFKXoZ5yKlEoSiEldDq8J3+FVsK7f9nx44m5JBmy8cFvYQMvfmsqVb5vWuZtpDTViLq4FfZOg1WDTfN3h1/Ok9g1GsG3L1eneUVPvt5ynubf7GT5kevqOYunlKoeqyiFXOTy5dz++BN+baWh8tB3eK3ya9nbcOUACFwLb2wBQzLcPQN3zsLdtFdShrIbjj6QEAFlm0PvJXka/6GgCL7YGMiJ61FU8HJgQruKNKvggRAiT4+jPJ4qM64ozygpJSFDhxG1dxf/G+zMr0P/xtYiG5MGxd6BGXUgOUNCsHICrwDwrASeAeBV2fS7jQvs+Ra2fQqvrYfSjfO8D5tO32bK5nMEhydQv4wb77WvSFVf5zw9jpI1lSgU5RlmCAvjQof2BFnHcuu7UQysOSR7GwbthhvH0pJCgOnM4VHf5PWJpsRi6wqDdoIm769cpxiMLDl0je+3XSQiPoVO1bwZ36aCmi2vgKhEoSjPuJitW7nx1kjWN7bmzZm7cLR0zPuDnFwOqwZB13lQ7eW833+a2CQ9c3dd4ae9VzAaoW/9Uoxo7o+LnWW+HVNRExcpyjPP8YUXoEML2u9N4s/VX+XPQap0h+LVYdtnpjOMfOJgbcG4NhXYOa453Wr68Ou+IBp8tZ2P15wmSD3VXeDUGYWiPENS4+I42qYJicYkIud+TNWSz+Pn5IdG5OF3wqA98NuL0PJj01DbAnDxTizzdl9hzYmb6I1GWlb0YmDj0tQt7apueuchdelJUYqIy9vWkDJ8Aj+/oGFLbQ0OFg5Udq9MVY+qVHWvynMez+Fq7fpkB1nS25QwRh4He4+8CTwb7sYmsejANRYduEpEfAqVvR0Z2Lg0HZ7zxlKnLpA8KZUoFKUIudKlCynCyLmv3+BU6ClOhp3kYuRFUqXpwTYfe5/0xNG0RFNKOJTI2QHCLsLMulC7P3SYmg89yFqSPpXVx2/w894gLt2Nw8vRin71/ehTtyTOtuo+Rm6pRKEoRUjEggXc+eJLSq9dg3X58gAk6BMIjAhMTxwnQ09yJ+EOdhZ2TGs+jXrF6+XsIBvGwZFfYNgBU2kQMzAaJbsuhvLL3iD2XAzDxkJL91q+jGtdASdbC7PEVJipRKEoRYghIoKLTZri2rcvXu++88h212KuMWrHKIJjgvmi0Re0K90u+weJD4Pvq4NfI3hlafpi/Z07SL0BS1+fJ+hBzp27HcPPe4JYffwGJVxt+bFfbfw97Qs0hsJOjXpSlCJE5+qKfbOmRK9bh9TrH9mupGNJfmv3G9U8qvHO7ndYdHZR9g9i5w6Nx8KFTabnMQCp13Pttde50qED0euyUdE2D1Us5sjXPaqxZHA9YhL1dJ25jx3n7hZoDM8qlSgU5Rnl3LUrqWFhxO3dm2U7R0tH5r4wl1YlWzH58GS+O/pd9st/1xsKTiXgrw/BaCRq9WpSgoOx8Pbm5vh3uDN5CtKQjbky8lAdP1fWvtWIEq62vPHbYebsuqzKmT8hsycKIURbIcR5IcQlIcSETNb3EUKcTHvtF0JUM0ecilLY2DdpgtbVlejVfz62rZXWim+afkPP8j355fQvfLjvQ/TGR5+JpLOwMQ2TvfUvxqOLCZs5C5vq1Smz5k9c+vQh4tdfuT54MIbIyCfvUA74ONuwcmh92lcpzlebzjFm2QmS9KpCbW6ZNVEIIbTATKAdEAD0FkIEPNAsCGgqpawKfA7MK9goFaVwEhYWOHV8kdgdO7L1Qa3VaPmw3ocMrz6ctZfXMnL7SBL0CY8/UNpDeJEzv8Bw5w4eY8YgLC0p9tGHFJ80iYTDRwju0ZOk8+fzoFfZZ2upY8YrNRjXujx/nrjJy3P/4Xa0ecqqF3bmPqN4HrgkpbwipUwBlgKdMzaQUu6XUv73r/wA4FvAMSpKoeXUrRvo9cRs2Jit9kIIhlQbwif1P2H/zf0M/GsgkUmPSTIaDakNPyD8WCp2AT7Y1X0+fZXzS90otWghMiWF4F69idm06Um6k2NCCEa0KMe8vrW4dDeOTjP2cvxawZ7dPAvMnSh8gOsZ3oekLXuUAUCm/9KEEIOFEEeEEEdCQwt+PmBFeRpZV6iAVUAlolevztF23ct357tm33Eh8gL9NvXjRtyNLNtHbD9PaooWj9IXHpqP26ZaNUr/sRLrihW5MWYsd6d+i0wt2MtArSsXY9WwhlhZaHh53gH+OBpSoMd/EokpqWafp8PciSKzZ+8z/YsIIZpjShTvZrZeSjlPSllbSlnbw6PgnhRVlKedc5euJJ05Q9L5CznarkXJFvzY+kfCk8Lpu7Ev5yMyv3RkiIggYv58HJo1wMYxDnZ++VAbnYcHpX6bj/PLLxP+449cHzKU1OjoTPaWfyoUc2Dt8EbUKunC2yv+ZdKGs2b/AM5MYkoqey6GMmXzObrO2keViVto/s1OTt8o2L9XRuZOFCFAxkdCfYGbDzYSQlQFfgI6SynDCyg2RXkmOHZ8ESwscnxWAVDDswYL2i5AIzS8vvl1ph6ZSmB44H2jiMLnzsWYlITHOx9A7Tfg6HwIfTipCEtLin86kWKffkr8gQME9ehJ8sWLT9K1HHOxs2TBgOd5rX4pftwTxBvzDxOTlI2b9vkoSZ/K/kthTP3rPD3m7Kfqp1vo+/Mh5u2+ggAGNiqNPtVIt9n7+f3gNbOM4DLrA3dCCB1wAWgJ3AAOA69IKc9kaFMS2A70k1Luz85+1QN3inK/kLfeIuHYccrt3IGwyPkTy7fjbzPp4CT2huzFIA2UdipNu9LtaGtTm+SXBuDYqSPekybdewjPqzL0WQHWmZc6Tzh2nJBRI5HxCfj9sRKr0qUffXCjEfZMhcQIaDgaHLxyHH9mlhy6xkd/nqaUmy0/vVaH0u52ebLfrOhTjVwNj+fS3TgCb8Vy4Eo4x69FkZJqRCPgOV9n6pdxo14ZV+r4uWJnpQMgIj6F0ctOsPtCKF1r+DCpaxVsLXV5Ht9T+2S2EKI9MA3QAr9IKScJIYYASCnnCCF+Al4CrqZtYsisIxmpRKEo94vdvp2QYcPxnTULhxbNc72fqKQotl7byqagTRy5fYTBGw00PQ1nZg6lxfMv42nrCadWwuo3waMivLIcnDK/7ai/cYNLrdvgNmAAnmPHZH5AfaJpvu7AtYAwDcetPxwavAXWTrnux38OXAln6KKjGCXM6lOThv7uT7xPMJ0lXA6N49Jd0+vinTguhcYRHBaPIe1ylxBQxduJemVcqV/WjTp+rjhYPzqJG42SGTsu8d3fFyjnac+sPjXx93TIk3j/89QmivygEoWi3E/q9Vxs1hzbmjXx/WF6nuzz5pnDRPV4jQMN3fiuURQCQZ1idWhXuh3NDFrc/3wLrBzg1ZWmM4xMXHtjACk3Qii7efPD5cLj7pqq1N44Cm0mQfm2sP1/cGYV2LiaSpzXGQgW1k/Uj+sRCQz47TCXQ+OZ2DGAvvX9crWfm1GJ/LD9EnsvhRISmch/H61ajaCUqy3+nvbpr3KeDpTxsEs/Y8iJvRfDGLX0OIn6VL7s9hydq+ddqRSVKBSliLvz1WQiFi+m3O5d6Fxcnnh/IaPHEL97N2X/3kqINoZNQZvYGLSR4JhgANwsnagYF0n55BQq1n6TChW7UsqxFDrNvQ/HyOXLuf3xJ5RevQrrSpXu7fzuOfi9h2kE1Us/QaUX7627edw0cdLl7eDoC83fg6q9QJv7SzGxSXpGLz3BtnN3ebVeST7pWBkLbfZu4UbEpzBrxyUWHLgKEloFeFLey4Fyng74e9rj526LlU6b69gyczs6ibeWHONwcCR965Xiwxcr5ckxVKJQlCIu6fx5gjp3weuDD3Dt++oT7SvxzBmCX+qO+7CheIwcmb5cSklgRCBH7xzlfMR5zoed4lLUFQxpJwtWWiv8nf2p4FqBCi4VqKT1wbrr8PsvP13ZCcv6mc4Uei8Fn5qZB3FlF/w9EW4eA/cKpifEK3Z49Jzfj5FqlHy95Txzdl2mfhk3ZvWpmeXUq3HJBn7ac4Wf9gSRkGLgpZq+jGpVDl+XgpnfW59q5Ost55m3+wpVfZ2Y+UrNJ55bXCUKRVEI6vYSAKVX/fFE+7k2cBBJp05R9u+taB2yvk6ujw/jyspXuRD6L+fLNuGcrQMXIi8QmWx68O3jZRLfWEuuzB1Lg+gw/LZ+jnAvD68sA+eSWQciJQSuM51hhF8E3zrQaqKpom0urToWwoQ/TlHc2ZqfX6v90H2AJH0qiw9eY+aOS0TEp9C2cjHGtSmf5/cLsmvLmduMW/EvGiH4tmc1WlbK/c1+lSgURSFi4SLuTJpE6TV/Yl2hQq72EX/oENf6vYbn+PG4DXgjexsZUmDdSPh3CdTsh2w/ldCUKE6HnebWkgXU/PUg7/TXElxMUBwdDcq0o0GJZtQtXhcnq2zctE41wInFsPMriL0JvZZAxfa56h/A0auRvLnwKMn6VKa/UoPmFTwxpBpZdfwG3/99kRtRiTTyd2d8mwpUK+Gc6+PklWvhCQxdfJQzN2OY378OzSp45mo/KlEoioIhMtI0T8Urr+D13kM1OB9LSsnV3q+gv3WLsls2o7HOwY1kKWHHF7B7Cvi3gh7zwcoBw91bXGzWAstK8RztVYd/3Hw5ePswsfpYNEJDFbcq1PeuT4BbAM5WzjhbOeNo5YiTpRMW2gdGCekTYV5zMCTB8IOgs8pxH/9zIyqRQb8d4dztGF5vUJpdF+5yOTSear5OvNO2Yp6NkHqc0FmzAPAYNizLdkn6VBYduMrrDfzQZfP+yoNUolAUBYCQkaNIOHKEcrt25viZitjtOwgZNoxin36Ky8s9cxfA0d9g/RjwCoAuc2DDWK4tvEiK0ZOyO/9BaDQYjAZOh51m/8397L+5n1NhpzBK40O7srOww8nSCSerey+PlCReO/A7xZp/DI1G5y7GNAkpBgYv+5MjkWsoZuzMhFaNaFPZ6+ERWvkkNTqai42bIPV6/FauwKZy5qPH8opKFIqiABC7YwchQ4fhO2smDi1aZHs7aTQS1KUrxuQkyq5fn6sH99Jd/BtWvAYpcaCzJtJhILdnr6L0qj+wDniwgDTEpMRwPfY60cnRxCTHEJUcRXRyNFHJUcSkxKT/Hp0czc24m9gaU5kUFkWTwQdz/YCelJJl55cx5fAU9EY9JRxKsLDdQtxs3HLf7xyKXLKE259+hsbWFqtKlSi1aGG+JqlHJYq8f7RPUZSnmn2jRmjd3YlevTpHiSJmw0aSL1zAe+o3T5YkAMq1gv4bTfcUGo3Fwd6f2/PWELNpc6aJwtHSkcpu2fs2HRwdzLjtbzGcVPpveJ23eq7FQpOzeOP18Xy6/1M2BW+isU9jepXqwtsHP2D4tuH80uYXbC0KZmRT1Oo/sapQAZdXXuH2J58Qu2kTju1zf+8lt8xd60lRlAJmmqeiI7E7dmKIiMjWNvq7dwn94QesKlXCsV0O5tXOSvFq0HsJlKiDzsUFu3r1iNmy5YlrGfk5+bG400petvXj1+QQ+q/tya24W9ne/kLkBXqt78WWq1sYVe0tPj5cAo/OY/je9U3ORZxj7M6x2ZvU6QklX7pE0smTOHfrinP3l7CqVIk7X3+DMTEx34/9IJUoFKUIcurSBQwGYtZveGQbY1IS0Rs2cG3QYC41a44+JATPt99GaPLnY8OxXVv0166RHBj4xPuy0lrxYcdFfB2VzKWoy3Rf150d13Y8drs1l9bQZ0Mf4vRx/FRvGq2nHyTytwWg0VBq2zk+qvcR+27uY+L+iflenC9q9WrQ6XB88UWEVkux99/DcOsW4T/9nK/HzYxKFIpSBFlXKI915cpE/Xl/RVkpJQnHjnHro4+52LgJN98eR/KlS7gNHkSZDeuxb9Qw32Kyb9kStFpiNm3Omx1aO9G20QesCLmBr86OkTtGMvnQZPSpD58NJBmS+GT/J3y470Oe83iOJQFTcBnxFfGHD1N80v9wefllYrdupYtXq/QZAL8/9n3exJkJaTAQvXYt9k2bonMz3ROxrVMHh3ZtCf/pJ/Q3Hyqyna9UolCUIsqpa1eSzwaSdO4c+hs3CJ01i8tt23L1lT5Er1+PQ4sWlJz/K/7b/sZz9OisK7zmgby4/JRy7RpXX+9PxMJFGJOTocarlHCvxMLrN3ilfE8WBS6i36Z+hMTem7joasxV+mzsw6qLqxj03CC+1fYmut8QjAkJlPrtN5xfegnnl7ohU1KI3rCBN6u+SY/yPfj59M8sDlycV92/T9zevaSGhuHcret9y73Gjwfg7jff5MtxH0WNelKUIsoQGcmlJk3ROjtjSJsV0rZuXZy6dMGx9Qto7PK/7PaDolau5NaHH+H3x8ocDwWVUnJ90GDi9+8HoxGdpydugwbhXNcXze+dodn7/F26Fh/v+xiAzxp+RqpM5ZP9n6DT6Piy4RdU2nCW0O+nYx0QgO/MGVgUK5a+/ytduyGEoPSqP0g1pjJ251h2XN/B102/po1fmzz9O4SMGk3C4cOZDmEO/WEGYTNnUmrRQmxrZ1lIO8ceNepJnVEoShGlc3HBuUcPNA4OeIwaif+2v02z0HXtYpYkAfcuP8Vu3pLjbWP/2kr83r14TXiXkvPnY1myJHcmTeLygI+JiKmHced3tHKuyPKOy/Fz8mPMzjGM2zWOss5lWd5yAaWnriZ02vc4duhAqcWL7ksSAM4vvUTS2bMkBQai1WiZ3GQy1T2r896e9zh8+3DOgr15AqbXhH+XPbTKEBlJ3PbtOHXsmOnoMreBA9AVK8btL74ouCllpZTP3KtWrVpSUZTC6eqAgfJiqxek0WjM9jap8fHyQtNm8nLnLtKo16cvjztwUAa/2leerVBRnq9WXoaPbSdTExNliiFFTj82XX5/9HsZfy1YXu7cRZ6tWEmG/fTTI49riIyUgc9Vlbf+Nyl9WVRSlOy0upOsv7i+PB9xPnvBGo1S/tpByk8cTa+/PpIy1ZC+OnzhInm2QkWZeO7cI3cRtX69PFuhooxYvjx7x8wm4IjM5DNVnVEoivJUcWzbBv316ySdPZvtbcJmz8Zw+zbFPv4Iobv3eJhd3ecptXABJX/7DasSxbmzIYhLLZoRu+h3hlcaxMDU+oT07I3+xg1KzJuL24ABj3ygTevsjEOrVsSsXYsxJQUAJysn5rSag43OhqFbh2ZvGO7lbRC8B1pPMs2nse97WNILkkxzYkevWoVVQKUsa3E5tm+PTc2ahH43jdTY2Gz/nXJLJQpFUZ4q9i1bgk5H7ObsjX5KvnyZ8F/n49S1K7Y1My9Jblf3eUqtXE/JFyVWdgnc+fIrLrVsxdX+b6B1dsZv2TLsGzd+7LGcXupGanQ0cdu2pS8rbl+c2S/MJtGQyJC/hxCdHP3oHRiNptLoziXh+cHQYSp0+NY0t8ZPrUg6+DdJZ8/i3LVblnEIIfB6/31SIyMJmzX7sXE/KZUoFEV5qqSPftr8+NFPUkpuf/4/NLa2eI57O+sdW9lj13cipRpepdQn/bCuWBGHFi3wW74MqzLZG9FlV78+Ou/iRP2x6r7l5V3K832L77kee53XN7/OmfAzme/gzCq4fQqafwi6tLku6gyAfmsgIZzoKUNAp8XxxQ6PjcWmSmWcXupGxMKFJF8JAiA5NTlb/cgplSgURXnqOLZrm63LT7GbNpFw4AAeo0elP2+Qped6gk9tbG/Mp+Ss7/Cd/v1j59PISGg0OHftRvy+fQ89y1CnWB1mtJxBTHIMfTb0YdrRafd/cBtSTFO5elWB53rcv2O/Rsj+W4m+YolD8Xh055bAY5IkgOfo0Wisrbkz+StO3D1B+z/acyr0VLb7k10qUSiK8tRxyMblp9S4eO58NRnrypVxefnl7O1Yo4F2kyHuNuz5NlexOXXtClIS9eefD61r4N2A1V1W09m/Mz+f/pnua7tz4u4J08pjv0FkELT8xBTHA+JOXSU1UeLUqDJseQ/WjABD1mcIOnd33IcNI37XbmbPGYSNhQ3e9t656ldWVKJQFOWpo3V2fuzlp7AZMzCEhlLsk48R2hzMF+1bG6q+DP/MhIigHMdm6euDbf16RK9ajTQ+XPrc0dKRTxt8ytxWc0lOTabfpn5MPjCJhF1ToGQDKPdCpvuNWrUKrbs79u+tgqYT4MQimP8ixN7JMp7Ers2566al11+JzGk2I1+q26pEoSjKUyn98tOZhy8/JV24QMTChTh3745N1ao533mriaDRwh8DIexSjjd3fqk7+pAQEg4demSbBj4NWN15NT0r9GTR+aW85KLlUO3emc7pbYiIIG7nLpw6dUJYWkLz96DHb3DnNPzYHE6thBtHIeYWGO89OxGeGM6QnSNY0tqGYmGp2KzZmeO+ZIdKFIqiPJXSLz9tuf/yk5SSO599jtbeHo+xY3K3c0dv6DwTwi7C7Aawc/JjL/PcF1urlmgcHR+6qf0gOws7Pqw6jF/CYhEWdgw4/jWf//M5cSlx97WLWb8eDAacunS+t7ByF3hjCwgt/DEAfmwB31aEzz3g2wDif2zOsBXtuBNznTfrt8Suuj9hP0zHcDM4B3+I7FGJQlGUp5LW2Rm7+vWJ2bT5vstPMevWkXDkCB5vj0Xn4pL7A1TpBiMOQ6UXYecXMLshBO3O1qYaa2ucXnyR2L/+IjUmJuvGe6ZSJy6aP1r/TL+Afqy4sIKua7uy98be9CZRq1ZjXaUK1uXL379t8aow4hAM3gW9l5qG0zYaQ0rpJoyyjOO8MZGp4THU2D8br+L70RKL/tRe8ppKFIqiPLUc27ZBHxKSfvkpNSaGO1O+xrpqVZy7d3/yAzh4Qfdf4NU/wKiH3zrC6iEQH/bYTZ1e6oZMTiZmw6NLtRN1HQ79CNV6Y1O8OuPrjGdh+4WmB/T+HsrI7SO5dGgryefO4fRAAcB0FjbgXR0qtIM6A0lt/j4TXGw5SBKfN/6CJmOD4L0QrN7/h7JL52LTtFPu/hZZUIlCUZSn1r3RT5sAU0G81PBwin38cd7Oi+HfCoYdgMbjTPcDZtSGYwtMD8g9gnVAAFYVK2Z9+WnnV6afzd5LX1TNoxorOq5gZI2RHLp9iDU/jCZVp0HfvO5jw5RSMungJLZe3cr42uPpWLajaYWVA3iUR5RrAdaO2epyTqhEoSjKUyv98tPmLSQFBhK5eDEuvXthUyVnlWWzxcIGWn4EQ/aCRyVY+xbM7wB3z2XaXAhhKhR4+jRJ588/3ODuOfj3d3h+EDiXuG+VldaKQVUHsaHjGlqes+BQOei4vTczT8wkXh//yBBnnpjJigsrGPjcQPpV7vdE3c0JlSgURXmqObZtiz4khJARb6F1dsZj1Kj8PaBnRXh9A3SaAaGBMKcR/P0pJD18L8LxxQ4ICwui/vjj4f1s/xws7aHR2EceyuLASaxik2n25mc08W3CnH/n0H5Ve5aeW/rQdKuLAxcz9+RcXir3EiNrjHzibuaEShSKojzVHFq2AJ0O/Y0beI4bh9bJKf8PqtFAzb4w4ojpKeq938K0KrDjC0i4N8+4zsUF+1YtiVm7Lr1QIADXD8G59dBgJNg9+rmG6FWr0Xl4ULpVZ75p+g2/t/+d0k6lmXRwEl3XdGXr1a1IKdlwZQNfHfqKliVb8mG9Dx9ZuDC/qEShKMpTTevsjGPr1tjWr3f/8NGCYOcOXWfD4J3g1xh2TYZpz8HWjyHuLmB6piI1Koq47dtN20hpKvxn5wn1hj5y14awMOJ278apS+f0irfPeTzHr21+ZUaLGeiEjrE7x9JrQy8+3PshdYrVYXKTyeg0ukfuM7+YPVEIIdoKIc4LIS4JISZksl4IIaanrT8phMi8PKSiKM8sn2+nUvKXX/L2BnZOeNeAXoth6D9Qvi3s/8GUMDa9i11ACXTFMxQKvPQ3XN0HTd8BK/tH7jJ63XpITTWVBMlACEHTEk1Z2Wklnzb4lLCEMMq7lmd68+lYaa3ys5ePZNapUIUQWuAC8AIQAhwGekspz2Zo0x54C2gP1AW+l1JmOTxATYWqKEq+Cr9suhz171JAEHq3LmHbgvD/eysWa1+GlFgYfvhehdgHSCkJ6tQZja0tfsuWZnkog9EAUCBnEk/rVKjPA5eklFeklCnAUuDBc8vOwIK0CZgOAM5CiOIFHaiiKEo6t7KmJ7tHHodar+FkdxSkJHpid7hzClp8lGmS0N+5Q+Ty5YQMHUbyxYsPnU1kRqfRmeVy030xmPXo4ANcz/A+BNNZw+Pa+ADZmEpKURQlHzmXhA5TsWw8DtuzLxN1/DZutaohKpsmHpJGI0lnzhC3YydxO3eml0238PbG9fXXceraxYzBZ5+5E0Vmt+4fvBaWnTYIIQYDgwFKliz55JEpiqJkl2NxnAe/w83x44krPQG5bRtxO3cSt2s3qWFhoNFgU6MGHm+PxaFZMyz9/Qt85NKTMHeiCAEyPoniC9zMRRuklPOAeWC6R5G3YSqKomTN4YVWaBwcCBllmmlP4+CAfeNG2Ddrhl3jxk9Wl8rMzJ0oDgPlhBClgRtAL+CVB9qsBUYIIZZiuiwVLaVUl50URXmqaKytKfbxxyQFBmLftCm2NWsgLCzMHVaeMGuikFIahBAjgC2AFvhFSnlGCDEkbf0cYCOmEU+XgASgv7niVRRFyYpTxxdx6viiucPIc+Y+o0BKuRFTMsi4bE6G3yUwvKDjUhRFUUzMPTxWURRFecqpRKEoiqJkSSUKRVEUJUsqUSiKoihZUolCURRFyZJKFIqiKEqWVKJQFEVRsmTWMuP5RQgRClw1dxy54A6EmTuIAlTU+guqz0VFYe1zKSmlx4MLn8lEUVgJIY5kVgv+WVXU+guqz0XFs9ZndelJURRFyZJKFIqiKEqWVKJ4uswzdwAFrKj1F1Sfi4pnqs/qHoWiKIqSJXVGoSiKomRJJQpFURQlSypRFDAhRFshxHkhxCUhxIRHtGkmhDghhDgjhNhV0DHmtcf1WQjhJIRYJ4T4N63PhXpyKiHEL0KIu0KI049YL4QQ09P+HieFEDULOsa8lo0+90nr60khxH4hRLWCjjGvPa7PGdrVEUKkCiG6F1RseU0ligIkhNACM4F2QADQWwgR8EAbZ2AW0ElKWRnoUdBx5qXs9BnTxFRnpZTVgGbAVCGEZYEGmrfmA22zWN8OKJf2GgzMLoCY8tt8su5zENBUSlkV+Jxn42bvfLLu83///idjmsWz0FKJomA9D1ySUl6RUqYAS4HOD7R5BVglpbwGIKW8W8Ax5rXs9FkCDkIIAdgDEYChYMPMO1LK3Zj68CidgQXS5ADgLIQoXjDR5Y/H9VlKuV9KGZn29gDgWyCB5aNs/HcGeAv4AyjU/x+rRFGwfIDrGd6HpC3LqDzgIoTYKYQ4KoToV2DR5Y/s9HkGUAm4CZwCRkkpjQUTnllk52/yLBsAbDJ3EPlNCOEDdAXmPK7t087sc2YXMSKTZQ+OT9YBtYCWgA3wjxDigJTyQn4Hl0+y0+c2wAmgBVAW2CqE2COljMnn2MwlO3+TZ5IQojmmRNHI3LEUgGnAu1LKVNPJcuGlEkXBCgFKZHjvi+lb9INtwqSU8UC8EGI3UA0orIkiO33uD3wlTQ/1XBJCBAEVgUMFE2KBy87f5JkjhKgK/AS0k1KGmzueAlAbWJqWJNyB9kIIg5TyT7NGlQvq0lPBOgyUE0KUTrtZ2wtY+0CbNUBjIYROCGEL1AUCCzjOvJSdPl/DdAaFEMILqABcKdAoC9ZaoF/a6Kd6QLSU8pa5g8pPQoiSwCqgbyE+O84RKWVpKaWflNIPWAkMK4xJAtQZRYGSUhqEECMwjYDQAr9IKc8IIYakrZ8jpQwUQmwGTgJG4CcpZZbD755m2ekzplEw84UQpzBdlnlXSlkYSzQDIIRYgmn0lrsQIgT4BLCA9P5uBNoDl4AETGdUhVo2+vwx4AbMSvuGbSjs1VWz0ednhirhoSiKomRJXXpSFEVRsqQShaIoipIllSgURVGULKlEoSiKomRJJQpFURQlSypRKIWSEKKnEOJ1c8dhDkKI1kKI0eaOQyk6VKJQCquewOvmDsJMWgOjzR2EUnSoRKEoTwEhhE1RPLZSOKhEoRQ6Qoj5wEtAUyGETHtNTFvXWQhxRAiRJIS4LYSYIoSwyLDtRCFEmBCiblq7RCHE3rQSI55CiD+FEHFCiEAhRIsHjhsshPhGCPFR2r7jhBCLhRBOD7RzFULMFULcSYtjvxCi7gNtpBBirBBimhAiFFPVXIQQHYQQW9MmxIkRQhwQQrTOGD/wNlAqQ9/np63bKYRY+cBxmqW1qZL23i/tfR8hxAIhRBSwLrtxK0WTKuGhFEafAyUBZ2BY2rIQIURPYAkwF3gfUyXaLzF9IRqXYXtbTBPnTAHigenAQiAZU/nrWcA7wAohRAkpZUKGbXtjKr0xCCieto+fSJtgSghhBfydFtt4TPMQDAX+FkKUk1LezrCv8cBuoC/3vrSVxvTB/Q2mEi7tgE1CiCZSyn1pxyqHqdJu17RtQrP7h8vgG0y1l3oAqTmMWylqpJTqpV6F7oWpyNrODO8FcBX49YF2bwCJgFva+4mYSno3zdBmWNqyjzMsC0hb1i7DsmBME9XYZ1jWB9MHeqW09wOAFKBchjY64DLwdYZlEjj+mD5q0rbdgqlG1n/LvwGCM2m/E1j5wLJmaceqkvbeL+396gfaZStu9SqaL3XpSXlWlMd0lrE8rfKuTgihA7YD1kCVDG1TgD0Z3l9K+7k9k2UPTii0VUoZl+H9KkxJqk7a+1bAUSAoQwwAuzCVnc5ow4OdEEL4CiF+E0LcwDTLnx7TzevymfT5STx47JzErRQx6tKT8qxwT/u58RHrM87/ECvvn0EvJe1n1H8LpJQpaVVOrR/Yz31TWkopE4UQcZguQ/0XRz1MH/APuvzA+zsZ3wghNJhKkDtgqrZ6CdOlsc8Az8w69QTuPPA+J3ErRYxKFMqz4r+5iwcDxzNZH5RHx7nvAzttxJA98N98EhHAEUzX9x+U/MD7B0s3+wM1MF3u2vzAMbIjCbB8YJnrI9o+eOycxK0UMSpRKIVVCvd/2z8P3AD8pJQ/5uNxXxBC2Ge4/NQN04fukbT32zBdKrompbyb2Q6y8F9CSP9gFkKUAhpimp/kPw/2/T8hQJMH483msZ8kbuUZpxKFUlidAzoLIbpg+oC8iWnY6EIhhCOm0UspQBmgC9Bd3j96KbcSgQ1CiK8xXW76GtON4bNp6xcAQ4CdQohvMM3U5wY8D9yWUn73mD6FAFOFEB9hugT1KaYE+GA7r7Qn009jmjo3GFgNDBBCfIfpHkRzTPORZ8eTxK0841SiUAqrWZgu0/wCuACfSiknCiFiMA2NfQNIxfSBt5579yGe1FIgFvgZ0yWntWS4XCOlTBJCNMd0X+FTwAvTfY1DPDwF7H2klMlCiG7ATEyjukKASZhGLmW8Gb8cUxKYAngAvwGvSyk3CCHexzSKayCmaXVHp/3M0pPErTz71Ax3ipJNQohgTMNPxz2uraI8S9TwWEVRFCVLKlEoiqIoWVKXnhRFUZQsqTMKRVEUJUsqUSiKoihZUolCURRFyZJKFIqiKEqWVKJQFEVRsvR/2e6SCz5/ZWIAAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "# 可視化\n", "for n,u4_beta in zip(n_list,u4_list_n):\n", " plt.plot(temp_list, np.array(u4_beta), label='n={}'.format(n))\n", "\n", "plt.legend()\n", "plt.ylabel('$U_4$', fontsize=15)\n", "plt.xlabel('temperature', fontsize=15)\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "統計が十分でないため、データにばらつきがあります。しかしおおよそ相転移点である温度が1付近で4つのシステムサイズのデータが1点で交わっていることとがわかります。 \n", "Binder cumulant による相転移点の推定は数値解析の最前線でよく用いられる手法です。\n", "\n", "> 当然ですが、学術的な研究では十分な統計を取ることはもちろん、誤差評価(エラーバーの計算)なども真面目に行う必要があるでしょう。今回の計算は概要の説明だけに留めてあるため、正確な誤差評価等は省略しています。\n", "\n", "## 結言\n", "\n", "アニーリングを用いてMonteCarlo samplingを行う方法をご紹介しました。それを応用して統計物理学における相転移の計算例を示しました。ここからもわかるように、OpenJijはアイデア次第で様々な応用が可能です。" ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [] } ], "metadata": { "kernelspec": { "display_name": "Python 3", "language": "python", "name": "python3" }, "language_info": { "codemirror_mode": { "name": "ipython", "version": 3 }, "file_extension": ".py", "mimetype": "text/x-python", "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", "version": "3.8.5" } }, "nbformat": 4, "nbformat_minor": 4 }