{ "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": 5, "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": 6, "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": "iVBORw0KGgoAAAANSUhEUgAAAYoAAAELCAYAAADHksFtAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAABau0lEQVR4nO3dd3gUxRvA8e/cXS6XXgmEJCShV6X3rnQFpSMgKogi2BB7w15A4IcI0kREikpHBUR6h9Ck1wSSACG9t7ub3x8XkZKEEC65kMznee5Jbm92992IeTM7s+8IKSWKoiiKkheNrQNQFEVRSjaVKBRFUZR8qUShKIqi5EslCkVRFCVfKlEoiqIo+dLZOoCi4O3tLYOCgmwdhqIoyn3lwIEDMVLKcrduL5WJIigoiJCQEFuHoSiKcl8RQlzMbbu69aQoiqLkSyUKRVEUJV8qUSiKoij5KpVjFIqiKP/Kzs4mIiKCjIwMW4dSYhgMBvz9/bGzsytQe5UoFEUp1SIiInBxcSEoKAghhK3DsTkpJbGxsURERBAcHFygfWx660kI8YMQ4poQ4lgenwshxFQhxDkhxD9CiIbFHaOiKPe3jIwMvLy8VJLIIYTAy8vrrnpYth6j+BHoms/n3YBqOa+RwIxiiElRlFJGJYmb3e3Pw6a3nqSU24QQQfk06QX8JC210PcIIdyFEL5SyitFEc+Oya+QdeUqWo0dWq0OrdYOrVaPVqO3fK+xQ6vRoBU6tBodwk4HdnZgZ/le6OwQ9nrQ2SGub7fDzuCAwcEVg4MLWnt7hF5ved34vfqHrChKCVXSxyj8gPAb3kfkbLstUQghRmLpdVCpUqVCnSz5j7+oFGFZnyO/rpY552VNJq3ArBOY9DrQ65AGe4S9PRoHR3QOjuicXLBzcELv5ILOwQmNwR6ht0cYDJbv7Q0Iez0ag8GyX85XyzEcLO8NhutfhcbWnUlFUQrq0qVL1K5dm/HjxzNu3DgA1q1bx8svv4zJZGLEiBG89dZbRXb+kp4oCkxKOQuYBdC4ceNCrcYUP2Y8p6IukpqZQlZ2KllZaWQb0xBkodVkoRFZaEU2IuelNYPWbMZOarCTEp2U6BFoJejMEp0EnRm00oQwG5GmbLKzjRiNJkzZJswmickEZrNAmgSYBRiz0Jiysc9OR28EfQro4yV6I9hnc9NXvbHwPy+TToPJXodZr8Ost8Nsb4eoUA6vJq3wadYGh3r10Dg4FP4EiqJYzdixY+nWrdv19yaTidGjR7Nhwwb8/f1p0qQJPXv2pHbt2kVy/pKeKCKBgBve++dsKxJPPNb/tm1SStKyTMSlZpGQlk18WpbllZpFcoaRpIxs4jOM179PyjCSnJ7zNSObTON/fQ8XBx11q7rxgL8b9fzdeNDPDX8XEFlpkJUMmSmQlUp2RgJJ6TEkpceSnJFAUmYCyVlJRGelkGRMJcmYRkZWKuaMJGRWJhhBGAXCDEh7NMKAEAY0Uo8w68AEIlsiss2ILBOaTCOaLCPaLCO6rGy02RnYZ0l8z8fgtO8kl76bg1kryKxcEbfGzSjXtDUODRpgV6FCUf3oFaVUCwsLo1u3brRu3Zpdu3bh5+fHqlWrcCjAH2MrV64kODgYJyen69v27dtH1apVqVy5MgADBw5k1apVZTZRrAbGCCGWAM2AxKIan8iLEAInex1O9joCPO9+/4xsE2GxqfwTkcjRiET+iUhg3s4wskyWBOLuaEc9Pzce9HenTsVyVHALwLucPV7Oerz0BfjPk5EIcaEQH3rz17hQSIoEcjpXupyXsx6cy4OzT85XyyvbyYsLnpU4FBvO5b2bkf+cIvBSJFV/W07W4uUAZHu7YWhQH6+GzbArXx6th8dNL41ef/c/IEUpRh+tOc6Jy0lWPWbtiq58+GidO7Y7e/YsixcvZvbs2fTv359ly5Zx5coVFi5ceFvbtm3bMnXqVFJSUvjqq6/YsGEDEydOvP55ZGQkAQH//Q3t7+/P3r17rXNBubBpohBCLAbaA95CiAjgQ8AOQEr5PfAn0B04B6QBT9sm0sIz2GmpWcGVmhVc6d/Y8h82y2jmTFQyRyIScpJHIjO2nsdkvvmOmYOd1pIwnO3xdtL/972zPZXLOdE0yBMngxtUrG953So7AxIjIOUqpERBcpTl67+v+IsQvg/SYrADaggNNeoPhqcnYXb24XzCeUIi9xJ2YDOZh//BPyyJGnu2YtqwNddr1Tg53ZA43NF5eKCrWBFDjRrYV6+BPrASQqu17g9YUe4TwcHB1K9fH4BGjRoRFhbGe++9x+uvv57nPuPHj+fVV1/F2dm5mKLMna1nPQ26w+cSGF1M4RQbvU5DXT836vq5WfpJWHoe566lEJ2cSUxKJrGpWcSmZBKbkkVMahZXEjM4djmR2JQsjDkJRacR1A9wp2UVL1pU8aZhoDv2uht+EdsZwLuq5ZUfUzYkX4E938O+WXBsOZpWL1Ot5Riq1R0CdYcgn5SEJoYSEhXCqpN/c/L8bpzTJE0MNWjt9ACVKYdMTMQYH48pPgFTXDxZ586TffUqmC29J2Fvj33VqtjXqIGhRnXsq1fHvkYNdJ6F6KopSiEU5C//omJvb3/9e61WS3p6OhMmTMi3R7F3716WLl3KG2+8QUJCAhqNBoPBQKNGjQgP/2+eT0REBH5+fkUWu7D8Li5dGjduLEtrmXEpJYnp2RyLTGLn+Rh2nY/laEQCZgn2Og1NgjxpWdWLllW8qVvRFZ32Lmc3xZ6HjR/BiVXg4gsd34MHB4Hm5p7A1dSrLDu7jOVnlnMt/Ro+jj70qdaH3tV6U8Hpv7EMc0YGmefPk3nmLJmnT5N55jQZp89gio293kZbzhuHOnUp//Zb6AMD7+nnoyi3OnnyJLVq1bJpDGFhYTzyyCMcO2Z5tnjixImkpKQwfvz4Ah9j/PjxODs7M27cOIxGI9WrV2fjxo34+fnRpEkTFi1aRJ06BU+Euf1chBAHpJSNb21b0scolFsIIXB31NO6mjetq3kDkJSRzd4Lcew6H8Ouc7F8ve40cBoXg45WVbx5o2sNKpcrYNfVqwr0/wku7YH178Kq0ZaeRpdPoXL7680qOFVgdP3RPPfAc2yL2MavZ37l+yPfM/OfmbT1a0u/Gv1oVbEVWoMBhzp1cLjlH7AxJobMM2fIOH2GzDNnSNm0ibD+A/D7dipOTZta6aelKKWTTqdj2rRpdOnSBZPJxDPPPHNXSeJuqR5FKRSdnMnuC7HsPh/Dn0evkm0y83GvuvRp6Hd3D/ZJCceXw9/jIeESVOsCnT4Gn5q5No9IjmDZ2WWsOLuC2IxYfJ18aVmxJd4O3ng7eOPl4IWXwev69446x+vxZF26RPjzo8gKD8d3/Hjc+/S2wk9CUUpGj6IkupsehUoUpdzVxAxeXnKIvaFx9KpfkU8fq4uLoWAVI6/LzoB9M2HbN5CVAvWfgMCW4BEE7oGWW1Q3PMCXbcpmU/gmlp1Zxun408RnxCO5/d+Zg84BT4MnXg5eBLgE8Er1Z8l6+3NSd+3Ca8Rwyo0dqx4MVO6ZShS5U4lCJYqbmMyS6ZvPMWXjWfzcHZg6qAH1A9zv/kCpsbD1Kwj5AczZ/23X2oN7JUviuPVVrgZGIUjITCA2PZaY9BhiM2Jv+j4mPYZ/ov+hglMF5nachXHyTBIWL8H54Yfw+/prNI6O1vgxKGWUShS5U4lCJYpchYTF8fKSw0QlZTCuSw1GtqmMRlOIGlPGLEgMh/iwXF4XITPxv7a+9eHJVeDgnn9sV0N4YeMLVHSqyNwucxFL1xL1xRfY16xBwPTp6mE/pdBUosidShQqUeQpMS2bt5b/w9pjV2lTzZtv+j+Ij4vBuidJj7ckjcgDsPYtyzMeQ1eAvUu+u+29spfRG0cT6BrI3M5z0e09QuTY19A4OuI/fToO9epaN06lTFCJInd3kyjUDeAyxs3RjumDG/JF73rsD4uj25TtbD59zboncfCAig2gyQjoNw8iD8KiAZCVlu9uzXybMbXjVMISwxi5YSSm5vUJXLQIYWfHxaFDSVr/l3XjVBSlQFSiKIOEEAxqWok1Y1rj7WzP0/P28+nvJ8g0mqx/slqPQu9ZcGk3LHnCMjCej5YVWzKlwxTOJZzj+Q3Pkx3sS9Cvv2CoUYPIl18mZuYsSmMvWFFKMpUoyrBq5V1YNaYVQ5sHMmdHKE/M3ktievadd7xb9fpCz2lwYTP8NswyxpGPNv5tmNR+EqfiT/H838+T4Wqg0k/zce3Rg+jJk7l2Q80bRSnNsrOzGTZsGPXq1aNWrVp88cUX1z9bt24dNWrUoGrVqnz55ZdFGodKFGWcwU7LJ4/V5dtBDfgnIoGBs/YQk5Jp/RM1GAw9voEz62DZcDDlXyO9fUB7JradyPGY47yw8QUyNCYqTpyAW98+xM37kcyzZ60fo6KUML/99huZmZkcPXqUAwcOMHPmTMLCwq6XGV+7di0nTpxg8eLFnDhxosjiUIlCAeDRBysyZ1gTQmNS6P/9biIT0q1/kiYjoMvncHI1rBwF5vxvdT0U+BBftf2KI9FHGL1xNBmmDHxeew2NkxNREyZYPz5FKSJhYWHUqlWLZ599ljp16tC5c2fS0+/8/5gQgtTUVIxGI+np6ej1elxdXW8qM67X66+XGS8qqoSHcl276uVYMLwZz8zbT78Zu/h5RLOCl/4oqBajITsdNn1iKVr4yP9ueljvVl2CumA0G3lnxzu8uOlFpnWchveoUVz7+mtSduzEuXUr68anlG5r34KrR617zAr1oNudb/0Upsx43759WbVqFb6+vqSlpTF58mQ8PT3LVplxpeRpEuTJ4pHNefKHffSfuZsFw5tRy9fVuidpO86SLLZPBJ0Bun0N+ZQW6VG5ByZp4r0d7/HK5leYMmgCdosXc+3rr3FqsVyVLlfuC4UpM75v3z60Wi2XL18mPj6eNm3a8PDDDxdTxP9RiUK5TV0/N359rgVD5uxlwMzdzHu6KY0CPax7ko7vgTEDdk+zJItOH+ebLHpW6YnJbOKDXR8wbtfbfDL2Za6+Oo6E5cvx6NfPurEppVcB/vIvKoUpM75o0SK6du2KnZ0dPj4+tGrVipCQEAICAoq1zLgao1ByVdXHmd+eb4Gnk56hc/ey42yMdU8gBHT+FBoPh11TYcsXliKE+Xi82uO81+w9tkZs5QuHLRgaNiD6f1MxpaRaNzZFKSavv/46hw8fvu01depUACpVqsSmTZsASE1NZc+ePdSsWZMmTZpw9uxZQkNDycrKYsmSJfTs2bPI4lSJQslTgKcjvz7fggAPR575cT9/Hb9q3RMIAd0nQv0hlhpSi/pbSoDkY0DNAbza6FXWXlzHmm5emGJiiJ07x7pxKUoJMXr0aFJSUqhTpw5NmjTh6aef5oEHHripzHitWrXo37+/KjN+t1QJD+tKSMti2Lz9HItMZELfB+jd0N+6JzCbYM8M2Pw5SDO0f8sy6K3Nu8rt1INTmX10NpO3VML/YCRV1q3FztfXunEppYIq4ZE7VcJDsSp3Rz0LRzSjWbAnY389wg87QolJybxtje+7ZTJLriVncOxKCgf9B8PovVClI/z9IcxsB5fynsXxYoMXGVRzEJ81uIjJbCR6ypR7ikVRlLypHoVSYBnZJsYsOsjfJy21oTQCPJ3s8XbW4+3831evf793scdkklxLziQqKYNryZlEJ2cQlZTJteQMYlKybko2859pSrvq5eDk77D2DUiKhEZPwcPjLfWjbmGWZt7f+T5Oc1fw+G5J0G+/qcKBym1UjyJ3ailUpUgY7LTMGNKIzaeucSUxg5iUzJxXFjEpmVy8lEpMchbp2bk/SOflpMfH1YCPiz21fF3wcTHg42qPj4s9H64+zpztFyyJotYjULkdbP4C9s6AU39YHtSr1++mmVEaoeGjlh/xTlISiUf+5sT412m49M+7W8VPUZQ7UolCuSt2Wg2d6+S/NkRqppHYlCyiUzLRagTlXe3xdrbHTpv3nc7z0alMWH+aU1eTqFnB1VKSvOvn8OAAWPMKLH8WDi+EHpMs63rn0Gl0fNrlG2bt7cNDv5xj+6KJtB2c97x0RVHunhqjUKzOyV5HJS9HGgV6UD/AHV83h3yTBMDgZpVwsNPyw47Qmz/wfRBG/G2ZHRV5EKa3gJ3/u2kqrV6r5+m3FhJd3gHTtHlsC91UFJelKGWWShRKieDuqKdPIz9WHrpMdPItRQk1Wmj6LIzeB9U6wYYP4NcnITP5ehMnB1dqfvAFFeIl6ya/yv6r+4v5ChSl9FKJQikxnmkVTJbJzM978niWwtUXBvwMnT+zjFvM7gjRZ65/XK5jZ/TNm9Bnu5E3fx/NsZhjxRS5ohSN2NhYOnTogLOzM2PGjLm+PS0tjR49elCzZk3q1KnDW2+9df2zzMxMBgwYQNWqVWnWrBlhYWH3HIdKFEqJUbmcMw/V9OHnPRfJyGNAHCGg5Rh4ciWkxVmSxck1OR8J/N5+D8dM6Lcbnv/7ea6mWvkhQUUpRgaDgU8++YSJuazBMm7cOE6dOsWhQ4fYuXMna9euBWDu3Ll4eHhw7tw5Xn31Vd588817jkMlCqVEGd4mmNjULFYdjsy/YXBbeG4rlKsOvwyBvz8CswlDjeq49+lN+73pOEclM/OfmcUTuKLko7Blxp2cnGjdujUGw83r2js6OtKhQwcA9Ho9DRs2JCIiAoBVq1YxbNgwAPr27cvGjRvveVVINetJKVFaVPaitq8rc7aH0r9xQP5TXd384em1lmcudkyCy4eg7w+Ue+klEv/4k3H7vHjdYyXP1HmGANeAvI+jlBlf7fuKU3GnrHrMmp41ebPpnf9qL0yZ8YJISEhgzZo1vPzyywA3lSDX6XS4ubkRGxuLt7f3XVzVzVSPQilRhBAMbx3M2WspbCtIIUKdPTz6P3h0KlzcCbPaocuOxPvZEfgdCKdOOMw4MqPoA1eUO8itzPidigLeidFoZNCgQbz00ktUrly5yGK3eY9CCNEV+B+gBeZIKb+85fNKwHzAPafNW1LKP4s7TqX4PPpgRb5ad4q5O0ItD+AVRKNhUL4u/DoUfuiCZ6evia/oy5htWTxXcQ3D6w2ninuVOx9HKdUK8pd/USlMmfE7GTlyJNWqVeOVV165vs3Pz4/w8HD8/f0xGo0kJibi5eV1T7HbtEchhNAC3wHdgNrAICFE7VuavQf8KqVsAAwEphdvlEpx0+s0PNkikG1nojkTlXznHf7l3whGbgX/JmjWvkT5jhVwuxhLl2M6vjv8XdEFrCiFdC89ivfee4/ExESm3FLnrGfPnsyfPx+ApUuX0rFjx3uuVmDrW09NgXNSygtSyixgCdDrljYS+HeJNTfgcjHGp9jIE80CMdhpbn8A706cy8HQldBsFC7Zf+BQK5jB2wU7Tv/FydiTRRKrohSloKAgxo4dy48//oi/vz8nTpwgIiKCzz77jBMnTtCwYUPq16/PnDmWcvvDhw8nNjaWqlWrMmnSJL788t4Xa7JpUUAhRF+gq5RyRM77oUAzKeWYG9r4An8BHoAT8LCU8kAuxxoJjASoVKlSo4sX81/XQCn53l1xlN8ORLDrrY54O9vfeYcbmYwwqRbp2jqEzT7D+uZ6zg1pzbSHphVNsEqJpYoC5q60lRkfBPwopfQHugMLhBC3xS2lnCWlbCylbFyuXAHvaysl2jOtg8ky5vMAXn60OnhwAA4p23Hr2Z1O+7I5c2QLR6KPWD9QRSnlbJ0oIoEb5y3652y70XDgVwAp5W7AABR+npdy36hSzpmOd3oALz/1h4DZiE97H7QGB4Zv0fLtoW+tH6iilHK2ThT7gWpCiGAhhB7LYPXqW9pcAh4CEELUwpIooos1SsVmhrcOJiYli9WHCzE05VMT/BqhC11BuVHP88CZLDJ27mbflX3WD1RRSjGbJgoppREYA6wHTmKZ3XRcCPGxEOLflcJfA54VQhwBFgNPydK42pKSq5ZVvKhZwYW5O0IL93Rp/cFw7QQeD9VDFxDAM5s0fBfyv3t+UlVRyhJb9yiQUv4ppawupawipfwsZ9sHUsrVOd+fkFK2klI+KKWsL6X8y7YRK8VJCMGINpU5HZXMjnMFeADvVnX7gM6A5sSvVHj7LXyjjfj8dZgdkTusH6yilFI2TxSKciePPuiLt7M9c+92qiyAgzvUfASO/oZz6xY4tGjOgB0wZ8dk1atQlAJSiUIp8ex1Woa1CGTL6WjOFuABvEuxaczfFcYHq46RmmmEBoMhIxFx5k8qvP02jpmS+qtPsfHSxmKIXlEKL68y4wBZWVmMHDmS6tWrU7NmTZYtWwYUTZlxm5fwUJSCGNw8kGmbz/HDzlC+6P3ATZ9lGk3sD41n8+lrbD59jQvRqdc/C/RyYnjLduDqD4cWYhi6HPeBA+m8eDFT135Dh5Ed0Gq0xX05ilIg/5YZP3bsGMeO3by+ymeffYaPjw9nzpzBbDYTFxcH3FxmfMmSJbz55pv88ssv9xSH6lEo9wVPJz29G/qz/GAksSmZXElMZ/G+S4z8KYSGH29gyNy9LNhzEX8PRz58tDabx7WnUaAH83eFYUID9QfB+U2QGIHPiy8inRzpuPIia0NV2TCl6Fm7zDjADz/8wNtvvw2ARqO5Xh1WlRlXyrThrYNYvO8SXaZsJybFslyqn7sDjzf0o0MNH1pU8cJR/98/6adbBTFm0SE2n7rGw/WfgG0T4MhidG1fx/flV9B+9jnzFn9Dl3e6Yqexs9VlKcXo6uefk3nSumXG7WvVpMI779yxnTXLjCckJADw/vvvs2XLFqpUqcK0adMoX768KjOulG1VfVwY1iKQqj5OvN2tJn+92pYdb3bg08fq8VCt8jclCYAudSrg62Zg3q5Q8KwMga3g8CKQEs+BAzFWqkC336NYc3K5ja5IKUusWWbcaDQSERFBy5YtOXjwIC1atGDcuHFFFrvqUSj3lY961S1wWzuthqEtAvl63WlOX02mRv3BsOoFuLQbEdiSoPc/RvfsSPbOmkzWpMfRa/VFGLlSEhTkL/+iYs0y415eXjg6OtK7d28A+vXrx9y5c4FSWGZcUYraoCaVMNhp+HFXKNTuBXZOcMjyP6ZLmzZkN3+QzpsTWbF3no0jVcqiwvYohBA8+uijbNmyBYCNGzdSu7ZlhYbSWGZcUYqUh5Oexxv4sfxgJPFGPdR5HI6vgMwUAGp88AV6kyD52+9JN955cFFRiltuZcYBvvrqK8aPH88DDzzAggUL+Oabb4BSWGa8qDRu3FiGhITYOgylhDh9NZkuU7bxZteajAqOgnndoNd0y/MVwNEPxqL5bS1HJz7FwB62WwFNKRqqzHjuSluZcUW5JzUquNCyihcLdodh9GtmGdg+/N994dpjPyDLoEMzYyEpWSk2jFRRSiaVKJQy4elWwVxOzGD9iWuWQoEXd0LcBQC07u7onx5EvXPZrPnlcxtHqiglj0oUSpnQsaYPlTwdmbczFB4cBEJjmSqbo9Zzr5HsacBj7mri0+JsGKlSFErjLfZ7cbc/D5UolDJBqxEMaxlEyMV4jiY7Q+UOcHgxmC0LImns7XF76QUCr5pYP/s9G0erWJPBYCA2NlYlixxSSmJjY3N92jsvajBbKTOSMrJp8flGutSpwKQ652HpMzB0BVTpCIA0m9nRtSUkJFHzrw2Uc/ezccSKNWRnZxMREUFGRoatQykxDAYD/v7+2NndXJEgr8Fs9cCdUma4Guzo1ziAhXsv8lbnh/AxuFmeqchJFEKjwe/Nt8kc/RZbJ79J349+tnHEijXY2dkRHBxs6zDua+rWk1KmDGsZRLZJsujANajXD079DukJ1z+v/FAvLj/gS9DKA1yOPG27QBWlBFGJQilTgr2d6FCjHD/vuURWvUFgzIBjy25qU+OdT3HIgn1fvmGjKBWlZFGJQilznm4VTExKJr9Hlwef2jc9UwHgX78ll9pUpcrGM4Sd3GujKBWl5FCJQilz2lTzpqqPM/N2XUTWfwIiD8C1m0tPN3j3a0xaOPH5uzaKUlFKDpUolDJHCMFTLYM4GpnIEY8uoNFByNyb2lQIrEV4jwYE74/k7E61uJFStqlEoZRJvRv64WrQMftQCtR/AvbNhtNrb2rT6vUJJDoJwr/8VM3BV8o0lSiUMslRr2NQ00qsO3aVyy0/At8HYPlIiD5zvY2nlx9XB7TF92w8J1f9ZMNoFcW2VKJQyqyhLQKRUrIg5BoMWAhaPSx5AjISr7fpOPoLrnhpiJ88FWk02jBaRbEdlSiUMsvfw5HOtSuweN8l0h0rQv+fID7U0rMwmwFwdfIgYfijeEalcXTeFNsGrCg2ohKFUqY93SqIhLRsVhyKhKBW0PVLOLMOtvxXRbbLkPc5W8mOrFk/YUpJtWG0imIbKlEoZVrTYE8eDHDn0z9OsOt8DDQZAQ2GwLYJcGIVAE56J4wvDMYpOZuj335i44gVpfipRKGUaUIIZg9thL+HA0/P28/m09HQYxL4NYYVoyDKsuxkj0df4WBdB8TiNWRHXbNx1IpSvFSiUMo8H1cDS0a2oFp5Z0YuCOHPk3Ew4Gewd4YlgyAtDnutPU5jRqIxmjk56SNbh6woxUolCkUBPJ30LHq2OQ/4uzNm0UGWnTVZkkViJCwbDmYTXVoPY1tDe7S/bybr4kVbh6woxcbmiUII0VUIcVoIcU4I8VYebfoLIU4IIY4LIRbl1kZR7pWrwY4Fw5vSoooXr/12hAWR5aHHN3B+E/w9HgedA2lDepClkUROmmjrcBWl2Ng0UQghtMB3QDegNjBICFH7ljbVgLeBVlLKOsArxR2nUnY46nXMHdaEh2v58P7KY8xMaQ2Nh8OuqXB0KY82HcofTQQZ6/8m/fhxW4erKMXC1j2KpsA5KeUFKWUWsATodUubZ4HvpJTxAFJKNZKoFCmDnZYZQxrxyAO+fLH2FFPshiMrtYBVY6iZmcmZbrVJddQSPWmyrUNVlGJh60ThB4Tf8D4iZ9uNqgPVhRA7hRB7hBBdczuQEGKkECJECBESHR1dROEqZYWdVsP/BjagXyN/pmwOY7LHe0hHT1gymEdqPcrSFpLUnTtJ3bPH1qEqSpGzdaIoCB1QDWgPDAJmCyHcb20kpZwlpWwspWxcrly54o1QKZW0GsFXfR7gqZZBTN2byMxy70JiON0S4tjWxIFUT0eufTNJFQxUSj1bJ4pIIOCG9/45224UAayWUmZLKUOBM1gSh6IUOY1G8OGjtXmhfRW+PO7OBYe6OB/8iQ5VurC4lZmMo0dJ/muDrcNUlCJl60SxH6gmhAgWQuiBgcDqW9qsxNKbQAjhjeVW1IVijFEp44QQvNG1JuM6V2dyYntE3AX6OlRiQ+1sMgN8iJ4yRRUMVEo1myYKKaURGAOsB04Cv0opjwshPhZC9Mxpth6IFUKcADYDr0spY20TsVKWjelYDddGfYiS7vjtX02QR2VWd3IhKzSUhOXLbR2eohQZURrvrzZu3FiGhITYOgylFMo0mlg26SWeSPuZye3f4oewhfz2ezW01+Kosn4dGgcHW4eoKIUmhDggpWx863Zb33pSlPuKvU5Lx8FvkI0Oz92H0Wnt2NErGOO1a8QvXGjr8BSlSKhEoSh3qYJfIAlB3emdsRUP8wPMt9uPY9s2xMyajSkx8c4HUJT7jEoUilII5R56CReRTtUIMwmZCZwb2AxzcjKxc+bYOjRFsTqVKBSlMPwbIys24ENxGHO2O98nbcf1kUeI+2kB2VFRto5OUaxKJQpFKQwhEE2fwy/rEjUyq3Aq8QBX+vVEms3EfDfd1tEpilWpRKEohVXncXD0YqJzIkjBK4dX49q3HwnLlpF5IdTW0SmK1ahEoSiFZWeARk9ROWwLdVzqEa/ZyfSA1gh7PdH/+5+to1MUq1GJQlHuRePhgOBZOwMauySWXD5IZKfeJK9fT/rRo7aOTlGsQiUKRbkXbn5Q6xHanvgbT3sPKvof4VVZG+nqRszMmbaOTlGswqqJQgjhZc3jKcp9oelz2GUk0MulGsmaf3D10fCXXwNStm7DlJJi6+gU5Z5Zu0ex38rHU5SSL7AllK9L7/ATmKSJHi0j2ez7AGRn888va2wdnaLcs0InCiGEWQhhyvlqFkKYgKAbvleUskEIaPosQVdP0NitOtuv/sHbr/Yh1smDwwuWMn71cVIzVXVZ5f51Lz2KGcACoJyUUiOl1AIXb/heUcqOev3B4E7vtEzCk8PRulwkqPejNIk+w29bTtL1f9vYdS7G1lEqSqEUOlFIKUcDc4CVQohh/262SlSKcr/RO0LDoXQ6twsXnRPLzi7D65HuaE1GFtTKRKfR8MScvby9/CjJGdm2jlZR7so9jVFIKXcAD2FZfGgDoLdKVIpyP2oyAoPZTA+DLxvCNpBZIxBdRV/KhWznz5faMLJtZX7Zf4nOk7ex5fQ1W0erKAV2x0QhhJghhHheCNFKCOFy6+dSyiwp5XvAWODToghSUe4LHkFQoxt9Lh0jy5zFH6F/4Nq1Gym7dqFPT+Gd7rVYNqolzvY6npq3n3G/HSExTfUulJKvID2KR4DpwDYgQQgRKoRYJYT4RAjR+98psVLKo1LK74syWEUp8Zo+S83Ea9R38mf64emkt20A2dkkb9wEQINKHvz+UmvGdKjKikORPDx5K38dv2rjoBUlf3dMFFLKAMAT6AC8DGwAygOvAEuBK0KIX4QQ/kUYp6LcHyp3AO/qfB6XCsDYqO/Q+lUkad3a603sdVrGdanBqtGt8Ha2Z+SCA8zbqWpDKSVXgcYopJQJUsptUsppUsqRUsrmUkoXoBbwFtAC2CuEqFiUwSpKiScENB1JwOUjTKj9LOcTL3CojoHUXbsxJSTc1LSunxurRreiS53yfLTmBIv2XrJNzIpyB/c6mH1aSjkJqA8YgfFWiElR7m8PDgS9Cy2Pr+PlBi+xsOJFMBpJ3rjxtqZ6nYZvBzWkQ41yvLvyKMsPRtggYEXJn1WezJZSxgFTgW7WOJ6i3NfsXaDtODizlqez7anRvCtR7nBx5eJcm+t1GmYMaUTLKl6M++0Iv/9zuXjjVZQ7sGYJj4uAtxWPpyj3r5YvQmBrxLo3+aj2cE7X90IbcpxLl47n2txgp2X2k41pHOjJy0sOqwFupUQpyPTYP4UQXwohBgkh6ggh8nrquhkQbt3wFOU+pdHC49+D0OK45mW6PDUerYQFM8eQlp2W6y6Oeh0/PN2Een5ujFl0SD1roZQYBelRuADPAQuBf4AUIcRBIcSPQog3hBAjhBDTgZeARUUYq6LcX9wDoMdECN9LYOZBTBV9qHowivd2voeUuRcxcLbXMf+ZplQr78xzCw6osh9KiVCQ6bFtpJQeQBDQC/gEOAM0AT4DZgFPYynn8XmRRaoo96N6/aBuH8TWr/Bp14x6F2HPyb+Ye2xunru4OdixYHgzgrycGD4/hP1hccUYsKLcrsBjFFLKS1LK36WUn0spB0op6wAOgB/gLKUcLaXMKrJIFeV+JAT0+Aacy+Nq2oAwS56JqcPUg1PZEbkjz908nfT8PKIZvu4Gnp63n8PhCcUXs6Lc4l6nxxqllFeklKqsuKLkxcEDHv8eexGK3tuRh84ZqOZRjTe2vcGlpLyfnSjnYs+iEc3xdNLz5Ny9HL+cWIxBK8p/1FKoilIcgtsiWo3BxSeKjP37mVx/PALBy5tfznNwG6CCm4FFzzbDxWDH0Ln7OBOVXIxBK4qFShSKUlw6vo9rfT8wS1w272VCuwlcSLzAezvfw2jOe2Ejfw9HFo5ohp1W8MTsvZy7ppKFUrxsniiEEF2FEKeFEOeEEG/l066PEEIKIRoXZ3yKYjU6e+xH/oDe1UTS4pm09G3Bqw1fZcPFDTy59kkuJFzIc9cgbycWjmiOEDBw1l7OXVNrcSvFx6aJIueZjO+wPNFdGxgkhKidSzsXLAUJ9xZvhIpiXaJCHVzbtyDtYjLGTd8xrM4wJrSdQHhyOP3W9GPu0bl59i6q+jiz+NlmAAyavYfz0SpZKMXD1j2KpsA5KeWFnBlTS7BMwb3VJ8BXQEZxBqcoRcFl+DsgBcnzJyBiz9M1uCsreq2grX9bphyckm/voqqPC4ufbYaUkkGz9nBBJQulGNg6Ufhx89PcETnbrhNCNAQCpJR/5HcgIcRIIUSIECIkOjra+pEqipXYV6+OPrgSSZfsYfmzYMrG28GbSe0nFah3Ua28C4uebY7JLBk0ew+hMak2uAqlLLF1osiXEEIDTAJeu1NbKeUsKWVjKWXjcuXKFX1wilJIQghcuz9KWpQO4/nDsOVLMJsRQhS4d1G9vAsLn21GtsnSswhTyUIpQiKvUgLFcnIhWgDjpZRdct6/DSCl/CLnvRtwHvi3f10BiAN6SilD8jpu48aNZUhInh8ris1lnj3LhUd7Ur5nNTwdt4LQgqMnOHqBoxfSwZP1OhOfpZ8hTZp4oUI7hgV1R2fvAhqdpb1Gw4W4TN5Ydgydzo5vBjbEz9PF8pm9Czi42/oylfuMEOKAlPK2CUO2ThQ6LOVAHgIigf3AE1LKXEtsCiG2AOPySxKgEoVyf7jw6KNoXV0JfKUDJF2BtBhIi4W0uJyvscRkxPOZlzt/OzlSLyOTj2PiqJpdgHW2NTp4fgf41Cr6C1FKjbwShc4WwfxLSmkUQowB1gNa4Acp5XEhxMdAiJRytS3jU5Si5NK1KzHTviO70mTsyvvk2sbbbGZSejzrz6/hs2Mz6efgwNCK7Xje72EchQ6kCcxGwmNTmLHpNI528EKbSnhuex/2zICeU4v5qpTSyKY9iqKiehTK/SDz/Hku9HiE8u++i+fQIXdsH5cRx5QDU1hxbgXlHcvzRpM36BTYCSEEAMciE3li9h5cHexYV/k3nE8vh7EnLbe0FKUA8upRlOjBbEUpzeyrVMG+enWS1q0rUHtPgycft/qYBd0W4GHw4LWtr/H8388TlhgGWNbgXjiiOUnp2Yw61xSMGXDgx6K7AKXMUIlCUWzItVtX0g8cIHnT5gLvU9+nPot7LOatpm/xT/Q/9F7dm6kHp5JuTKeevxs/j2jGWVmJnea6JG2fQWamevxIuTcqUSiKDbkPGIB9rVpEvPACUV9PQBZkoBrQaXQMrjWYNY+voUtQF2Yfnc1jKx9j06VN1PNzY/2rbTkTNBjXrGtMmDJRlSlX7olKFIpiQzpPT4KWLMbjiUHE/fADF4cMJTsyssD7ezt480WbL5jXZR6Odo68vPllxmwaQ5LxKk8/9TzpzpXolbGa3tN38tW6U2RkqxUBlLunEoWi2JjG3p4KH3yA3+RJZJ47x4XefUjetOmujtG4QmN+ffRXxjUeR8jVEPqs7kN4SiQOrV+gnjzNK7WSmbHlPI9+u0P1LpS7phKFopQQrt26Ebx8GXo/PyJeGE3Ul18hswq+aKSdxo5hdYaxtOdSss3ZLDy1EOoPBr0LLzlt5Menm5CSaaT39J18uVb1LpSCU4lCUUoQfWAggUsW4zF4MHE//kjYkKFkRRT8VhRAgEsA3YK6seLsCpI1AhoMhuMraF/RzPpX29K/cQDfbz1Pj6nbOXQpvoiuRClNVKJQlBJGo9dT4f338JsyhawLFwjt3Zvkv/++q2MMrj2YNGMaK86ugKYjwWyE/XNxNdjxZZ8HmP9MU9KzTPSZsYtxvx1hf1gcpfGZKsU61AN3ilKCZYWHE/nqWDKOHcPjyaGUHzcOodcXaN9ha4cRlRbFH4//gXbJYIjYD68eBzsDAEkZ2Uxcf5qlByJIyzIR6OVIn4b+9G7oh7+HY1FellJClchaT0VFJQqlNDFnZXFtwkTiFyzArmJF3Pv3x71vH3Te3vnut+HiBsZuGcuU9lN4yKSFn3pBr+mWW1E3SM00su7YVZYeiGD3hVgAWlT2om8jf7rVq4Cj3qaVfpRipBKFotznUrZvJ3buD6Tt2QN2drh26oTHoIE4NG58vYzHjYxmIz2W96Cic0XmdfkBZrS0VJZ9fjvk0h4gPC6NFYciWXoggktxaTjptXSr50vfRv40DfJEo8l9P6V0UIlCUUqJzAsXiF+yhMQVKzEnJ2NfrSruAwbi1qsnWheXm9rOPz6fiSET+fWRX6kVthfWvAxP/QFBrfM9h5SSkIvxLA2J4I+jV0jJNBLo5chnj9WjdbX8ezLK/UslCkUpZczp6ST9+Sfxi5eQcewYwtERtx498Bg0EENty9LzSVlJPPzbw3QK7MRnTd+FybUhsBUMXFjg86RnmVh//CpTN53lQnQqT7cK4s2uNTHYae8pfiklmUbzPR9HsR5VFFBRShmNgwPuffoQvPQ3gn77FdeuXUlcvZrQ3n24OGQoppQUXPWuPFb1MdaGriXGlAaNnoLTf0L8xQKfx0Gv5bEGfvzxYhuGtQhk3s4wHv12B8ciEwsde0hYHP1n7qbWB+t4b+VREtMLVrpEsQ2VKBSlFHCoV4+Kn39GtW1b8Rn3GmkhISQsXQrA4FqDyTZn8+vpX6HJCEDAvll3fw69lo961WX+M01JTM/m8ek7mb7lHCZzwe9KnL6azIj5++n7/W7CYtPo9WBFFu29xEPfbGHloUg1RbeEUreeFKUUCntiMMboaKqsW4vQahm9cTTHYo6xoe8G9Mufg3MbYewJsHcu2AGlhIwEcPAAID41i3dXHuXPo1dpEuTBpP71CfDMe0ptRHwakzacYcWhSJztdTzfrgpPtwrCUa/jWGQi7648xpHwBFpW8eKTx+pSpVwB41KsSt16UpQyxHPoELLDw0nZtg2AIbWGEJcRx9rQtdB8FGQmwpHFdz6QlHBmPczuCBOqQoTlDzAPJz3fPdGQSf0f5NSVZLr9bzu/hYTf1iOITcnkozXH6ThxK7//c4WRbSqz7fUOjO5Q9fq027p+biwf1ZJPH6vL0chEuk3ZzqS/TqsSIyWI6lEoSikks7M593An7KtUodIPc5FS0nt1b3QaHb/2+AUx5yHITILR+0GTy9+L/yaIrV/C5UPgXgmy0sAzGIZvuGl6bUR8GmN/PcK+0Di61CnPF70fQK/TMHd7KLO3XyAty0i/RgG80qkavm4O+cZ9LTmDz/84ycrDlwn0cuTjXnVpV72ctX88Sh5Uj0JRyhBhZ4fHoEGk7tpF5rlzCCEYUmsIp+JOEXLtgKVXEXsOzm+8eUcp4dSfMKs9LB4AaXHQcxq8eBAeHm95uvvYspt28fdwZPGzzXmne002n4qm8+RttPt6M5P/PkPrqt789Wo7vur7wB2TBICPi4EpAxuwcEQztEIw7Id9jF50kKgktfiSLakehaKUUsa4OM6174Bbn974fvghGcYMOi3tREOfhvyv7QSYUg/K14Ghyy0J4vSfsOVLuPoPeARB29fhgQGgtbMc0GyyJJC0OBizH/S3j0mcvJLEG0v/wdVBx+tdalI/wL3Q8WcaTczceoFpm8+h12r4pv+DdKlTodDHU+5M9SgUpYzReXri2qMHiatWY0pKwqAz0K96PzaHbyY8PQqaDLf0KPbOhJltYMkTkJlsKfMxJgQaDPkvSQBotND1C0iKgN3Tcj1nLV9X/vdkRaY8UeWekgSAvU7LSw9VY8OrbfH3cODjNSfuaoaVYj0qUShKKeYxZDAyLY2E5csBGFhzIFqhZfGpxdDoadDqYe0bkJUKj32fkyAG35wgbhTUGmr1hB2TIenybR+bpZkR60fw+rbXrXYNgV5OvPJwNSIT0vn7ZJTVjqsUnEoUilKKOdSpg0PDhsQvXIQ0mfBx9KFzUGeWn11Oit4B+syB3nMsg9r1B4G2AAUAO31sKVu+8ePbPjoVd4ro9GgORB3gROwJq13Hw7XKU9HNwPxdYVY7plJwKlEoSimX21TZ1OxUVp1fBbV7wQP9CpYgrh8wGJq/YJleG3ngpo92RO4AwEHnwMKTBS8Tcic6rYYhLQLZdT6WM1HJVjuuUjAqUShKKefy8MPoypcnfsHPANQrV4/65eqz8ORCTOZCPqvQ5jVwKgfr3rYMhOfYEbmD2l61/ysbkh5jjUsAYGCTStjrNPxYQnsVZrNk6sazbDl9zdahWJ1KFIpSylmmyg60TJU9fx6AIbWHEJ4czraIbYU7qMEVOr4P4XvhuGX8IzEzkSPRR2hVsRVP1HyCbHM2v53+zVqXgaeTnl71K7LiYCSJaSWvNtRX604xacMZPlh1HHMpG3RXiUJRygD3/v0Rej3xCy23gx6q9BAVnCrw88mfC3/QBkOgfD3Y8CFkp7Pnyh7M0kwb/zYEuQXRxq8NS04vIcuUZaWrgGEtg0jPNvHbgXCrHdMa5u8KY+a2C9T1c+VSXBrbz1mvJ1USqEShKGXAv1NlE1auwpSUhE6j44maT7Dv6j5Ox50u3EE1Wuj6OSSGw+5p7IjcgYvehXre9QBLryUuI451Yeusdh11KrrRNMiT+bvDSsxU2XXHrjJ+zXEerlWepc+3xMtJz897Cl6d936gEoWilBG3TpXtXa03DjoHy1TZwgpuCzUfQW6fzM6IbbTwbYFOYxkYb+HbgipuVfj5xM9WrQo7rGUQ4XHpbD5l+7GAAxfjeXnJIR70d+fbQQ0w2Gnp1ziAjSejuJKYbuvwrEYlCkUpI26dKutm70anwE6sD1tPpimz8Afu/AlnNGaiM+Jo7fffynlCCAbXHszJuJMcvHbQCleQc7o65fF1MzB/d1ihjyGl5Iu9XzBx/0RSslIKdYwL0SmMmL8fXzcDc4c1xkFvWYDpiaaVkMDifcV/e6yoKm3YPFEIIboKIU4LIc4JId7K5fOxQogTQoh/hBAbhRCBtohTUUqDW6fK9qjcg5TsFLZHbL+Hg1Zme412ALTWud/00SOVH8HN3s2qU2XttBqGNA9k+9kYzl0r3FTZPVf2sOjUIuafmE/PlT1ZG7r2rn7JxqRk8tS8/Qgh+PHppng521//rJKXI22rlWPJvktkm8yFiq8w9oXG0WfGLsLj0qx+bJsmCiGEFvgO6AbUBgYJIWrf0uwQ0FhK+QCwFPi6eKNUlNLj1qmyTSs0xcvgxR8X/rin4+4w6KmRbabc5q9umi7roHOgT7U+bLy0kciUyHs6x40GNglAr9Mwf9fdjwVIKZl6cCq+Tr781O0nfBx9eGPbGzy74VlCE0PvuH9alpHhP+7nWnIGc4c1Jsjb6bY2Q5oHci05k43F+CT55A1nuBSXjvcNSctabN2jaAqck1JekFJmAUuAXjc2kFJullL+myL3AP7FHKOilBq3TpXVaXR0C+7G1oitJGUlFeqYyVnJHIk5RmvfZnBpN5xYedPng2oOQiBYcmqJFa7AwsvZnp4PVmTZwQiSMu5uquymS5s4FnuMUQ+OooFPAxZ2X8h7zd7jRMwJeq/uzdSDU0k35j6+YDSZeXHRIY5GJvLtoIY0qOSRa7uONX2o6GZg4d5Ld31thbH7fCy7L8Qyqn2V67fArMnWicIPuPFGXkTOtrwMB9bm9oEQYqQQIkQIERIdHW3FEBWldLl1qmyPyj3INmfz98W/C3W8vVf2YpRGWjd8HsrXhb8+gOz/yoJXcKrAw4EPs+zMMtKyrXdb5KmWQaRlmfgtJKLA+5jMJr499C1BrkF0C+pBptGEVqNlQM0BrH58Nd2DuzP76GweX/U4W8K33LSvlJIPVh9n46lrfNSrLp1ql8/zPFqNYGDTSmw/G0NoTGohr7BgpJRM/vsMPi72DG5WqUjOYetEUWBCiCFAY2BCbp9LKWdJKRtLKRuXK6cWOlGUvOg8PXHt3v36VNk6XnUIdA3kzwt/Fup4OyJ34GznzIPlG0CXzyHxkmUti+2T4OwGSLrCkJqDSc5OZvX51Va7jrp+bjQK9GDB7rACP+D2Z+ifnE88T+/gEbSbsI16H/7FY9/t5KM1x9l9JotRdd7lh84/YNAaeHHTi7y46cXrt8ymbznPor2XeL5dFYY2z2Wo1JgFC/vD7umA5faYViNYtLdop8ruPh/LvtA4XmhfBYOd9XsTAHdR4KVIRAIBN7z3z9l2EyHEw8C7QDsp5T1Mz1AUBcBjyBASV64kYflyvJ56iu7B3fn+yPdEpUZR3invv5RvJaVkR+QOmvs2x05jB5XbQZtx8M8vcGHL9XYPOnhSt4InC0Om0D/NiMa3HnjXADtDznrciZB8FZIvQ9IVy9fkqzd/7+gNlZpBQHPLV7cAnmoZxIuLD7H1TDQdavrkG2u2KZvvDn9HkEt1pqy2x14LT7UK4nB4Aov3XWLezjAAfFzseTDgLXw8drAr8hd6rXyMFl69+XO/I10frMfrnavnfoJNH8PZ9RB1DJo9j4+rgc61y/PbgQhe61yjSH6JSymZtOEMFVwNDGxaNL0JsH2i2A9UE0IEY0kQA4EnbmwghGgAzAS6SiltP3FaUUoBh7r/TZX1HDqUHpV7MOPIDNaFrWNYnWEFPs65hHNEpUUxym/Ufxsfet/ySo+HqBMQdRwRdZTB0Qd4m2R2bniNNukZILTg6gdpMZDbLSmDO7hWBJcK4FMbEiPg8GLYP8fyuUtFevg35ayTB1s2x9Gh2sC8y6MDy88uJzIlEs2Vp/G3S2NW/+r4BVcHjZZsk5lTV5I5FB7PwYvxHApP4OKJWgjdKzhUWMMW0yIcA2FnFrRc4kiQWxBBrkEEuwUT5BZEcEIUgbunYfCuATGnLaVNAlswpHkga49d5c+jV+jd0PrDqzvOxRByMZ5PetUpst4ElIAV7oQQ3YEpgBb4QUr5mRDiYyBESrlaCPE3UA+4krPLJSllz/yOqVa4U5Q7S1q7lshXx+I/YzouHTow6PdBmKSJXx/9tcDHmHdsHpMOTGJD3w1UcMp/9blsUzZdlnWhurM/3wc+DlHHIeGSpbigSwVw8bW8XHO+2uWydKrJCNeOw6W9EL7H8jXJMkZh1jmg8W8MFeqBMcPSS8l5pWck0t0xHf+sbH66epXrK377N4UnV+W6Wl9sSiaHLiVwKDye6LRrdG6gITojnLDEMMKSwghNDOVK6pXr7YUEX6fyVIwLp5xbJbyrdMLb4M3sLbG46734bkB7vB29cbFzQdyw5nhhSSnpPWMXUYkZbH69Pfa6e08Uea1wZ/NEURRUolCUO5PZ2Zzr3AWZmYnvp5+yskIEX+//mlW9VlHZvXKBjjFi/QhiM2JZ0WtFgdrPPDKTaYen3dU57iQ28gKfzJjHQN/LNNedhegzoHcCg9v111RjMrO11xgWXZ1RDzbByc3Lsprfpk+gWmcYsPDuSq3nSM9K5dLivoTGHCe0+bOEmVK4GraVmOxkou0dSTfdvta3vdYebwdvqrlX49PWn+Jm71ao695y+hpPzdvPZ4/XZXAz6zxelleisPWtJ0VRbETY2VFpzmwix71OxAsv0LpPL74NFvwR+gcvNnjxjvunZqdy4NoBhtYaWuBz9qvRj1n/zGLhyYW83+L9ewn/Oi+/ymjq9WH48avseechXAw3335aeugssw4NxtFYh+GjfsLJSf/fhwY3+GMs/P4y9JwGd/mXvsO+2dQI3UWNRyZD42dyAloJvw1DDl1FakBjLsRfpt+cv2hRTUeHOg5Ep0VzLf0aGy5u4L2d7zG1w9S77mFIKZm84Qx+7g70axRw5x3u0X0z60lRFOuzr1KFoF+W4DViOBnLVzNlvo6j21YU6CnlvVf2YjQbbyrbcSeeBk96VO7BmgtrSMxMvJfQbzKsZRCpWSaWHbh5quxvIeG8u2kaQpvOd93ewePGJAGWdcPbvQmHfoZNn97dSSNCLD2SWj0ty8r+q1pnsHNCHF+Os96ZB8pXp0e11uw/FkSfqoMZ12QcX7f9mtcavcaW8C0sOLHgrq938+lrHIlI5MWOVdHriv7XuEoUilLGafR6fMaNo9KPP+KMnpe+v8Kxb8YjjcZ899sZuRNHnSMNfBpc35Z95QpxP/1E+pEjee43uNZg0o3pLDu7zGrX8GCAOw0quTN/98XrU2UX7A7jjRW7MHjt5KGATjSp+EDuO7d/GxoOg+0TYe+sgp0wIxGWPmMZS+k59eaeiN4RanSDk6vBZHkYcEjzQFKzTKw89N+kzsG1BtMxoCOTD07maPTRAl+rpTdxlgBPB/o0Kp7nj1WiUBQFAKdmTQlasYK9tbXo5vzKxaFPkhWee2G7f6fFNvNthkhKIX7JL1wcMpRzHToS9fkXXH7zLaQ59zpHNTxr0KRCExafWozRnH8yuhtPtQwiNCaVbWejmbn1PO+vOk7V6vtAk83LjfK5lSYE9JgENbrD2jfg+Mr8TyQl/D7WMgurz1xwyOXp7Lp9LLO+cqYINwhwp5avKwv3XrreWxNC8HGrj/Fx8OH1ba8XuIf198lrHI1M5MWO1bDTFs+vcJUoFEW5ztXblxMvdmFuH1cyz50j9LHHSVix8rZbUReunCB4XwSDf7jI2TZtuTp+PMb4eMq9/BI+r48jKyzseuHB3AypNYSrqVfZeGmj1WLvVteXci72jPvtH75Ye4qH6+mJ02yhV5VeBLsF57+zVgd9f4CAZrD8WQjNp0ji4YVwbKmlJ1KpWe5tqj4E9m5wzFLSXQjBkOaVOHkliYOXEq43c7N3Y0K7CUSlRvHBzg/ueMvv37GJQC9HejfIr4iFdalEoSjKTXpU7sH66mnEfP8Ohtq1ufL220S+8irG6GiSN28m8rVxpHcbxMurzbiFJ+D11DCCV66g8u9r8B41Cs8nn0RXvjxxP87P8xzt/Nvh5+xn1aqyep2GIc0CiUnJpF8jfyoG7UAief7B5wt2ADsHGLQYPCvDkifgai63g2LOwp+vQ1AbaDM272Pp7KHWI3Dq9+vlTHrV98PZXsfCW57UfqDcA7zS6BU2hW9i0alF+Ya4/ngUJ64k8VLHauiKqTcBKlEoinKL1n6tcdG78HvaXir9OA+fca+RvGkTZ9u0JWLUC6Tu2MGxxl58/5w/NTZvxmfcOAw1a16fuSPs7PAYMpi0PXvIOHUq13NoNVoG1xrMoWuHOHztsNViH9W+Cj8905RRnVxZfX4V/Wv0p6JzxYIfwNEThiwDexf4uS/E3/BL3ZgJS58GnQF6z7Ks8JefOr0hMwnOW3pNzvY6HmtQkd//uUJ86s3Lwz5Z+0na+7dnYshEjsccz/VwZrNkyt9nqOztRK/6d3FNVqAShaIoN9Fr9XQO7MzGSxtJN2fiNWIEwb8swfOZZ/D/fgZ+m9fzdftE/Ft3Rmhy/xXi0b8/wsGBuPk/5XmePtX64GnwZNqhadaLXaehbfVyfH9kBnqtnhH1Rtz9Qdz8LcnCmA4/94bUWMv2DR9YehmPTbc8MX4nlduBgycc+2/QfnCzQLKMZpbeMjtLCMGnrT/F28GbcVvHkZx1+zob645f5dTVZF56qHh7E6AShaIouehRuQfpxvTrFVQNtWtT/o3XcWnfngNxR8g2Z9PKr1We+2vd3HB//DGSfv8dYx7VnB3tHBlRbwR7r+5l35V9Vov9dNxp1oatZUitIXg7eBfuID614IlfLQPWi/pZftnv/R6aPW+Z0VQQWjuo3RNOr4UsSwXZWr6uNAr0YNG+S7cVMnSzd2NC2wlcSb3Ch7s+vGm84t/eRJVyTjz6YPH2JkAlCkVRctGofCPKO5bnj9DbFzTaHrEdB50Djcvf9gDvTTyGDkVmZxO/OO81ufvX6I+Pgw/TDk+z2jKe3x76Fhe9C0/VfereDlSpuWWA+/Ihy1TYCvWg08d3d4y6fSx1rM6sv75pSPNKhMaksut87G3N6/vU56WGL7Hh4gaWnP5v/Y4/jl7hTFQKLz9cHa3m3st/3C2VKBRFuY1GaOge3J1dkbuIz4i/vv3fabFNKzRFr9XncwSwDw7GuUMH4hcvwZxxeykLsJSzGPnASA5dO8Suy7vuOe4DUQfYGrGVZ+o+g6ve9Z6PR80e8OhUywB3nx8sg9R3I7AVOJeH48uvb+pW1xcPR7vbBrX/9VSdp2jj14YJ+ydwIvYEppzeRDUfZ3rU872Xqyk0lSgURclVj8o9MEojf4X9dX3bpeRLRKRE5Hvb6Uaew4Zhio8ncc2aPNv0rtabik4V+fbQt/fUq8g0ZTJ+13gqOlXkiZpP3HmHgmo4FF46BOXyKC+eH40Waj8GZ/6CDMsKggY7Lf0aB7D22FUafbKBrlO28eQP+3jt1yN8te4U83ddpL3nSzjr3Hll02ss2n+G89GpvGKj3gSoRKEoSh6qe1SnqnvVm24/7YjcAVDgsh2OzZpiX7MmcfPn55kE7LR2PP/g8xyPPX7bqnJ3Y8bhGYQlhfFhyw9xtLu9GqzN1O0NpkzLWEWOUe2qMLZTdTrXqYC/hwOJaVnsOh/D7G0X+GjNCd749TyRZ/pwOfUyn+39mBoVnOlWN//qvEVJFQVUFCVXQgi6B3dn6qGpRKZE4ufsx/bI7QS5BhHgUrBCdEIIPIcN48rbb5O6YyfObXJPMI9WeZQ5R+cw7fA02gW0QyPu7m/Y4zHH+fH4j/Su1puWFVve1b5Fzr8puPpbBsQfHACAh5Oelx6qdltTs1kSn5ZFdEom15KasuxCBpvEj5TzWsDHezbjZu+Gu7077vbut33vZu+GTlM0v9JVolAUJU/dK1sSxdpQyyyikKsh9K3e966O4dqjO9cmfUPc/Pl5JgqdRseo+qN4e/vbbLi4gS5BXQp8/GxTNu/veh8vgxevNX7trmIrFhoN1HkM9s6EtDjLsxp5NhV4Odvj5WxPzQrQutqrzJl/kb3XzrMt7QLxmfH5lj1xsXNhWc9l+DpbdyxDJQpFUfLk5+xHA58G/HHhD2p41CDTlHlX1WLBUnTQ84kniP7fVDLPncO+atVc23UL6sacf+Yw/fB0Hq70MNo7PdCWY87ROZyNP8u0jtOsM4BdFOr2gd3TLE9qN3yywLul7wuhzdcbae/pSfAvS9BVrEiaMY2EzAQSMhNIzEj87/u0aBKijhR6fYv8qDEKRVHy1T24O+cSzvHj8R+x19rfcVpsbtwHDkTY2+f7AJ5Wo+WF+i9wIfECf4b+WaDjnok/w6x/ZtE9uDvtAtrddVzFpmID8Ai6XvupIIwxMUSOew07f39kVhbhzz+POSUFJzsn/Jz9qONVh5Z+Leke3I0nTPaM2jqTt/ctwzE+99lU90IlCkVR8tUlqAs6oWPf1X00rtAYg85w18fQeXjg1rMniatWYYyLy7Pdw4EPU9OzJjOOzCDbnJ3vMY1mIx/s/ABXe1feavrWXcdUrISw9CpCt0JK7g8g3kiaTESOex1zUjL+307Ff+pUMsMuEvnyy8jsG34uUcdh/qOW5zwcPeGZ9ZaHBa1MJQpFUfLlYfCgpZ9lgLiNX5tCH8dz2JPIrCzilyzJs41GaBhdfzThyeGsOW+ZUmtOTeXqZ5+TsGLlTW1/OvETx2OP806zd/Aw5FLqu6Sp0xukGU6uumPTmOkzSNuzhwrvv4ehRg2cmjfD9+OPSd21myvjxyPT4uDPN+D7NhB1zFImfeRWy0OCRUCNUSiKckd9q/Vlz+U9tPMv/O0d+6pVcWrThvhFi/EaMQKNPvcH9tr5t6Oedz2+P/I9nYw1uDb2dbJCQ0Gjwc63Ak7NmxOaGMp3h77j4UoP0zmwc6FjKlbl64B3DcvtpyZ516BK3bWLmOnTcevVC7c+fa5vd3/8MbIvXSRmxvfow37Fu3qMZfnVDu/mO0BuDapHoSjKHXWo1IGdg3bi73JvK6p5DhuGKSaGpD/yHoMQQjCm/hhq7I7k0oCBmFKS8Z8xHX3lYCJfeZWM8HA+2PkBBp2Bd5u/e9frTduMEJZnKi7ugqTLuTbJjrpG5OtvoK9SmQoffnDztUWE4G1YjmtgGtEHdCRW+Rx6fFPkSQJUolAUpYAKMzZxK6dWLbGvVpW4H3/M8wE8c1oagd+u4YU/zJz10+D32xJcOnQgYNo0pMnE8eee5PjlQ7zV9K3CF/2zlTq9AZnrKnrSaOTya69hTkvDf8oUNI45Dw0mR8HKF2DOQ4iUKHw//xLHxo258tUM0g4cKJawVaJQFKXYCCHwePJJMk+fJm3v3ts+zzx/ntD+/UlatYqsYY/xQX8Ty2It6znog4LQf/wGhgtXeX+rNz2CexR3+PeuXHUoX++m2k//iv52GmkhIfiO/xD74CBLIcFfn4QpdeGfX6HVKzBmP5pGg/Cf9i12fn5EvDCarLCwIg9bJQpFUYqV26OPovX0vG0FvMTVqwnt2w9TXDwBc2bz4Ntf0MyvBXOPzSUtOw0pJZ9q1rKqvYGa+6NI+Nl6q+MVq7q9IWL/TYsipWzbRuzMmbj1eAg3wz6YVAsW9YewHZbxjNF7odNHlgWVAK27OwGzZoJGw6XnnsMYH5/X2axCJQpFUYqVxmDAY+BAUrZsITM0FHNGBlfef5/Lb7yJQ506BK9YgXMrS9HBMQ3GEJcRx6JTi1h6din7ru6j+itv49yxI1FffUXqPuutY3EvsqOukbRuPVFfTyBuwc/5/+Ku87jl6/EVln3Pn+Dy2Jex99ZSwWGB5QnugKYwcDG8dhq6fgFeVW47jL5SJfy/+w7jlatEjB6DOTOzKC4NAGGtGvAlSePGjWVISIitw1AUJQ/GmBjOdeiIc/t2ZF0KJ/P0abxGjqTcSy8idDdPxhy9cTSHrx3GJE3U9a7L7E6zMaekENZ/AKbERIKXLcXOt/jKb0ujkYzTp0k/dJj0Q4dIP3SI7Ms5g9N2dpCdjbCzw6XTw7j16YNTixa3rwQ4uyOkxyO9anFxxl4yE3QEDfbCvuMwqNsXnLwKHE/SunVEvvIqrt27U3HihDxXHSwIIcQBKeVtT1SqRKEoik1cfvsdElesQOvuTsWvv8K5bdtc252IPcGA3wfgoHNgec/l12deZV64QFi//uiDgwlc+DMa+7tcK6KApMlE6s6dpB08aEkO//yDTE8HQOfjg0PDhjg2qI9Dw4YYatQg88IFEpYuI3HNGsyJidhVrIhb7964934cu4o5q9Pt+R7WvUnU8QrEHdVQ8f1XcBv8XKFjjJk9m+hvJuH13HP4vPpKoY+jEoWiKCVK9uXLxM77Ea9nnr5jj2DBiQUEuATQPqD9TduT//6biDEv4vb44/h+/pnVp8rK7Gwix71O8vr1oNViqFkThwYNcGhQH8cGDdD5+uZ5TnNmJsl//03ismWk7toNQuDUsiXu/fri3K4tqb8vJOL9SbgPHIDv+PH3FqeUXP3gQxKWLiV45UoMNQqxdgYqUSiKUkpFT51KzPQZlP/gfTyfsN6CRTI7m8ixr5G8YQPlxo7Fc8jg/6as3qWsiAgSly8nYfkKjFevonV3RxqN2FUKIGjxYqv0hmR2NmkHDuLUvFmhj1FiE4UQoivwP0ALzJFSfnnL5/bAT0AjIBYYIKUMy++YKlEoStkhzWbCR40idecuAuf/iGOjRvd+zKwsIl97jeQNf1P+7bfwHDbMCpHm3MbatYuEpcvIPH2agJnfow8MtMqxraFEJgohhBY4A3QCIoD9wCAp5Ykb2rwAPCClfF4IMRB4XEo5IL/jqkShKGWLKSmJ0H79MKemWQa3y5cv9LFkVhYRr44lZeNGyr/zDp5PDrVipCVbXonC1tNjmwLnpJQXpJRZwBKg1y1tegH/TrheCjwk7ptn9hVFKQ5aV1cCpk3DnJZGxIsvkX31aqGOY87KIuLlVyxJ4v33ylSSyI+tE4UfEH7D+4icbbm2kVIagUTgtrljQoiRQogQIURIdPSdy/gqilK62FerRsWvviTz5EnOd+nKtUmTMSUnF3h/c1YWkS++RMrmzZbxjsGDizDa+4utE4XVSClnSSkbSykblytXztbhKIpiA66dOlF57VpcOncmdtYsznfqTNxPC5BZWfnuZ87MJOLFF0nZupUK4z+06qB4aWDrRBEJ3LhKu3/OtlzbCCF0gBuWQW1FUZTb6P398JvwNUFLl2JfsyZRn3/O+R6PkLR2ba6FCM2ZmUSMeZHUrduo8NFHeAwcaIOoSzZbJ4r9QDUhRLAQQg8MBFbf0mY18O+Ug77AJmnrqVqKopR4DnXrUGneDwTMnoXGwYHIV8cS1n/ATWU/zBkZRLwwmtTt26nwycd4DOhvw4hLLpsmipwxhzHAeuAk8KuU8rgQ4mMhRM+cZnMBLyHEOWAsUMLXPFQUpaQQQuDcpg3BK5bj+/nnGK9d49KTwwgf9QLpx45bksSuXfh+9ike/frZOtwSy+bPURQFNT1WUZTcmNPTiVvwM7GzZmFOSQEh8P3sM9x7P27r0EqEvKbHqqVQFUUpMzQODniPfBb3fn2Jm/cjhtq1ce3axdZhlXgqUSiKUuboPDzwGfuqrcO4b9h6MFtRFEUp4VSiUBRFUfKlEoWiKIqSL5UoFEVRlHypRKEoiqLkSyUKRVEUJV8qUSiKoij5UolCURRFyVepLOEhhIgGLto6jrvkDcTYOohipq65bFDXfP8IlFLetk5DqUwU9yMhREhuNVZKM3XNZYO65vufuvWkKIqi5EslCkVRFCVfKlGUHLNsHYANqGsuG9Q13+fUGIWiKIqSL9WjUBRFUfKlEoWiKIqSL5UoipkQoqsQ4rQQ4pwQItf1v4UQ/YUQJ4QQx4UQi4o7Rmu70zULISoJITYLIQ4JIf4RQnS3RZzWIoT4QQhxTQhxLI/PhRBias7P4x8hRMPijtHaCnDNg3Ou9agQYpcQ4sHijtHa7nTNN7RrIoQwCiH6FldsVielVK9iegFa4DxQGdADR4Dat7SpBhwCPHLe+9g67mK45lnAqJzvawNhto77Hq+5LdAQOJbH592BtYAAmgN7bR1zMVxzyxv+TXcrC9ec00YLbAL+BPraOubCvlSPong1Bc5JKS9IKbOAJUCvW9o8C3wnpYwHkFJeK+YYra0g1ywB15zv3YDLxRif1UkptwFx+TTpBfwkLfYA7kII3+KJrmjc6ZqllLv+/TcN7AH8iyWwIlSA/84ALwLLgPv6/2OVKIqXHxB+w/uInG03qg5UF0LsFELsEUJ0LbboikZBrnk8MEQIEYHlL68Xiyc0mynIz6Q0G46lR1WqCSH8gMeBGbaO5V6pRFHy6LDcfmoPDAJmCyHcbRlQMRgE/Cil9MdyW2aBEEL92yyFhBAdsCSKN20dSzGYArwppTTbOpB7pbN1AGVMJBBww3v/nG03isBy/zYbCBVCnMGSOPYXT4hWV5BrHg50BZBS7hZCGLAUVbuvu+v5KMjPpNQRQjwAzAG6SSljbR1PMWgMLBFCgOXfc3chhFFKudKmURWC+quteO0HqgkhgoUQemAgsPqWNiux9CYQQnhjuRV1oRhjtLaCXPMl4CEAIUQtwABEF2uUxWs18GTO7KfmQKKU8oqtgypKQohKwHJgqJTyjK3jKQ5SymApZZCUMghYCrxwPyYJUD2KYiWlNAohxgDrscyG+EFKeVwI8TEQIqVcnfNZZyHECcAEvH4///VVwGt+DcsttlexDGw/JXOmjNyPhBCLsSR775xxlw8BOwAp5fdYxmG6A+eANOBp20RqPQW45g8AL2B6zl/YRnmfV1ctwDWXGqqEh6IoipIvdetJURRFyZdKFIqiKEq+VKJQFEVR8qUShaIoipIvlSgURVGUfKlEodyXcirsPmXrOGxBCNFZCPGKreNQyg6VKJT7VX/gKVsHYSOdgVdsHYRSdqhEoSglgBDCoSyeW7k/qESh3HeEED8CfYB2QgiZ8xqf81kvIUSIECJDCHFVCPG1EMLuhn3HCyFihBDNctqlCyF25JQY8RFCrBRCpAghTgohOt5y3jAhxEQhxPs5x04RQiwUQrjd0s5TCDFLCBGVE8cuIUSzW9pIIcRYIcQUIUQ0cDRnew8hxIacBXGScioId74xfixPsgfecO0/5ny2RQix9JbztM9pUzfnfVDO+8FCiJ+EEAnAmoLGrZRNqoSHcj/6BKgEuAMv5GyLEEL0BxYDM4F3gCrAF1j+IBp3w/6OWBZL+hpIBaYCC4BMLOWvpwNvAL8JIQKklGk37DsIS+mNZwHfnGPMAfoBCCHsgb9zYnsdS2HDUcDfQohqUsqrNxzrdWAbMJT//mgLxvKLeyJgxrLIz1ohRFsp5c6cc1UDOmIpYQ2Fq4s1EUvtpX6A6S7jVsoaW6+cpF7qVZgXliJrW254L4CLwLxb2j0DpANeOe/HY6kn1e6GNi/kbPvghm21c7Z1u2FbGJaFapxv2DYYyy/0WjnvhwNZQLUb2uiwrPI34YZtEjh4h2vU5Oy7HkuNrH+3TySXVQCBLcDSW7a1zzlX3Zz3QTnvV9zSrkBxq1fZfKlbT0ppUR1LL+NXIYTu3xeWZSgNQN0b2mYB2294fy7n66Zctt26oNAGKWXKDe9XYElSTXLePwwcwFIi/t8YALZiKTt9oz9vvQghhL8QYr4QIhIwAtlYBq+r53LN9+KPW97fTdxKGaNuPSmlhXfO19t++ea4cf2HZHnzYjJZOV8T/t0gpczKqXJquOU4N62RIaVME0KkYLkN9W8czbH8gr/V+VveR934JmexptWAC5Zqq+ew3Br7GPDJ7aLuQdQt7+8mbqWMUYlCKS3+Xbt4JHAol89DrXSem35hCyEcAWfg3/Uk4oAQLPf3b5V5y/tbSzdXBRpgud217oZzFHRWUgagv2WbRx5tbz333cStlDEqUSj3qyxu/mv/NJZV4oKklLOL8LydhBDON9x+ehzLL92QnPcbsdwquiSlvNsV+v5NCNd/MQshAoFWwD83tLv12v8VAbS9ZVvnXNrl5l7iVko5lSiU+9UpoJcQ4jEsvyAvY5k2ukAI4Ypl9lIWUBl4DOgrb569VFjpwB9CiAlYbjdNwDIwfCLn85+A54EtQoiJWFYn9AKaAlellJPvcE0RwDdCiPex3IL6iNuXST0FlM95Mv0YECOlDMMyXjJcCDEZyxhEB3KWmC2Ae4lbKeVUolDuV9Ox3Kb5AcvtlY+klOOFEElYpsY+g2WFwAvA7/w3DnGvlgDJwFwst5xWc8PtGillhhCiA5ZxhY+A8ljGNfZx+xKwN5FSZgohegPfYZnVFQF8hmXm0o2D8b9iSQJfA+WA+VhWBfxDCPEOlllcI4BVwMs5X/N1L3ErpZ9a4U5RCkgIEYZl+um4O7VVlNJETY9VFEVR8qUShaIoipIvdetJURRFyZfqUSiKoij5UolCURRFyZdKFIqiKEq+VKJQFEVR8qUShaIoipKv/wOoy6rdVPW3eQAAAABJRU5ErkJggg==\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.9.3" } }, "nbformat": 4, "nbformat_minor": 4 }