{ "cells": [ { "cell_type": "markdown", "metadata": {}, "source": [ "# A3 - Classical ising model" ] }, { "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/en/A003-LargeScaleMC.ipynb)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In this section, we explain how to use core interface (core Python interface) of OpenJij and demonstrate simple calculation.\n", "\n", "Core interface is lower layer API than tutorials in previous section. It is assumed that the target reader has done though the previous OpenJij tutorials and is familiar with terms such as ising models and MonteCarlo methods.\n", "\n", "We envision a readership with the following objectives.\n", "\n", "* want to use OpenJij not only for optimization problem but also for more specialized application such as sampling and reserach activity.\n", "* want to set annealing schedules and directly improve the algorithms used." ] }, { "cell_type": "code", "execution_count": 1, "metadata": {}, "outputs": [], "source": [ "import cxxjij.graph as G\n", "# set problem size N = 100\n", "N = 100\n", "\n", "graph = G.Dense(N)\n", "# for sparse, use the following\n", "#graph = G.Sparse(N)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We set $J_{ij}, h_i$. In this time, we set the value generated from the Gauss distribution with a mean of 0 and a standard deviation of 1." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!pip install numpy # use numpy for random number generation" ] }, { "cell_type": "code", "execution_count": 3, "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "mu, sigma = 0, 1\n", "\n", "for i in range(N):\n", " for j in range(N):\n", " # to avoid a large Jij value, normalize in 1/N\n", " graph[i,j] = 0 if i == j else np.random.normal()/N\n", "\n", "for i in range(N):\n", " graph[i] = np.random.normal()/N" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "For longitudinal magnetic fields, both `graph[i]` and `graph[i,i]` can be used to access it. In addition, by definition of ising model, $J_{ij}$ and $J_{ji}$ are automatically the same value. Let us try output as follows." ] }, { "cell_type": "code", "execution_count": 4, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "0.5\n", "0.5\n", "-0.6\n", "-0.6\n" ] } ], "source": [ "graph[20] = 0.5\n", "print(graph[20,20])\n", "print(graph[20])\n", "graph[12,34] = -0.6\n", "print(graph[12,34])\n", "print(graph[34,12])" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## Set System\n", "\n", "Next, we define system for calculation. We can choose classcal ising model or transverse magnetic field ising model or other models.\n", "\n", "In order to create system of classical ising model, we choose `system.make_classical_ising`." ] }, { "cell_type": "code", "execution_count": 5, "metadata": {}, "outputs": [], "source": [ "import cxxjij.system as S\n", "\n", "mysystem = S.make_classical_ising(graph.gen_spin(), graph)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The first argument is randomly generated spins, and the second is the Graph. We can make a system of classical ising model with an initial spin configuration of `graph.gen_spin`. \n", "\n", "We can access system directly and can read values." ] }, { "cell_type": "code", "execution_count": 6, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[ 1. -1. 1. 1. -1. 1. 1. 1. -1. -1. -1. -1. 1. -1. 1. 1. -1. -1.\n", " -1. -1. -1. 1. -1. -1. 1. -1. 1. 1. -1. 1. 1. -1. 1. -1. -1. -1.\n", " -1. 1. -1. 1. -1. -1. -1. 1. 1. 1. 1. -1. 1. -1. -1. 1. 1. 1.\n", " -1. 1. -1. 1. -1. -1. 1. 1. 1. 1. 1. -1. 1. 1. 1. -1. 1. 1.\n", " -1. -1. 1. -1. -1. 1. -1. -1. -1. -1. 1. 1. 1. 1. -1. -1. -1. 1.\n", " 1. 1. -1. -1. 1. -1. -1. -1. 1. -1. 1.]\n" ] } ], "source": [ "print(mysystem.spin)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Various other systems including classcal ising models are available, which can be used for different purpose. A method of initialization differs slightly depending on the system to be used. We will introduce this later.\n", "\n", "## Run algorithms -Updater, Algorithm-\n", "\n", "After definition of System, we choose Updater and run Algorithm.\n", "\n", "### Updater\n", "\n", "There are certain updaters that can be used for the System. There are two main Updater for classical ising model, \n", "\n", "- [SingleSpinFlip](https://github.com/OpenJij/OpenJij/blob/ec41aecfbac7e4c895e1e7a1718f06eb7ffae0ba/src/updater/single_spin_flip.hpp#L40) (Update one spin at a time using the Metropilis Hastings method)\n", "- [SwendsenWang](https://github.com/OpenJij/OpenJij/blob/ec41aecfbac7e4c895e1e7a1718f06eb7ffae0ba/src/updater/swendsen_wang.hpp#L45) (Cluster update scheme using SwendsenWang method)\n", "\n", "### Algorithm\n", "\n", "Algorithm needs a **schedule list** to be executed. First, we begin with making schedule list.\n", "\n", "#### Schedule list\n", "\n", "Schedule list is a list of `(parameters, the number of MonteCarlo steps)`. The value you enter in the parameter part depends on the system. For example, in the case of classical ising model, parameter is the inverse of the temperature, $\\beta$. Let's take a look at the settings as follows." ] }, { "cell_type": "code", "execution_count": 7, "metadata": {}, "outputs": [], "source": [ "schedule_list = [(0.01, 10),(10, 80),(0.1, 30)]" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In above setting, 10 MonteCarlo steps at an inverse temperature $\\beta=0.01$, 80 steps at $\\beta=10$, 30 steps at $\\beta=0.1$, for total of 120 MonteCarlo steps.\n", "Since the inverse temperature is often increased by the geometric series, it is more convenient to use `make_classical_schedule_list` in `utility` as follows." ] }, { "cell_type": "code", "execution_count": 8, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[((beta: 0.100000) mcs: 20), ((beta: 0.199474) mcs: 20), ((beta: 0.397897) mcs: 20), ((beta: 0.793701) mcs: 20), ((beta: 1.583223) mcs: 20), ((beta: 3.158114) mcs: 20), ((beta: 6.299605) mcs: 20), ((beta: 12.566053) mcs: 20), ((beta: 25.065966) mcs: 20), ((beta: 50.000000) mcs: 20)]\n" ] } ], "source": [ "import cxxjij.utility as U\n", "schedule_list = U.make_classical_schedule_list(0.1, 50, 20, 10)\n", "print(schedule_list)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "In above example, we set calculations are performed in 20 steps each, varying the temperature in 10 steps from $\\beta=0.1$ to $\\beta=50$.\n", "\n", "#### Run Algorithm\n", "\n", "Next, we run Algorithm. By writing `Algorithm_[Updater]_run`, the calculation is performed in the specified Updater. The following example shows SingleSpinFlip." ] }, { "cell_type": "code", "execution_count": 9, "metadata": {}, "outputs": [], "source": [ "import cxxjij.algorithm as A\n", "A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "The process took an instant but a total of 200 MonteCarlo steps were computed during this time.\n", "\n", "> `A.Algorithm_SingleSpinFlip_run(mysystem, seed, schedule_list)`, we can keep the seed fixed. It can be used when we want our results to be reproducible." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can use callback to get the system at every single MonteCarlo step during Algorithm execution. In the case of classical ising model, we just need to create a funciton with a system and a parameter (inverse temperature) as arguments.\n", "\n", "The following example shows callback which record energy values." ] }, { "cell_type": "code", "execution_count": 10, "metadata": {}, "outputs": [], "source": [ "energies = []\n", "\n", "def callback_log_energy(system, beta):\n", " # graph is defined in Graphe module before\n", " energies.append(graph.calc_energy(system.spin))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We execute same Algorithm with this callback." ] }, { "cell_type": "code", "execution_count": 11, "metadata": {}, "outputs": [], "source": [ "# we take long schedule (a total of 20000 MonteCarlo step)\n", "schedule_list = U.make_classical_schedule_list(0.1, 50, 200, 100)\n", "A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list, callback_log_energy)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "Let us make the energy of the recorded system with MonteCarlo steps on x-axis, and energy on y-axis, we get following figure." ] }, { "cell_type": "code", "execution_count": null, "metadata": {}, "outputs": [], "source": [ "!pip install matplotlib" ] }, { "cell_type": "code", "execution_count": 12, "metadata": {}, "outputs": [ { "data": { "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEGCAYAAACD7ClEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA0N0lEQVR4nO3dd5wU9f348df7Gp2jHYj0JkXpCBbACioYMWqIioklUWMviYbYNfkqsUZ/xoItmhh7iyLSBASkSO9wlBM44DiOckc5rn1+f8zssnu3/XZvdm/fz8djH8zOzs68d2+Z98ynijEGpZRSySnF6QCUUko5R5OAUkolMU0CSimVxDQJKKVUEtMkoJRSSSzN6QDC0aJFC9OxY0enw1BKqYSyZMmSvcaYLF+vJVQS6NixI4sXL3Y6DKWUSigi8rO/17Q4SCmlkpgmAaWUSmKaBJRSKolpElBKqSSmSUAppZKYJgGllEpimgSUUiqJaRKoAZNX7aLg0DGnw1BKqSo0CcTYvsMl3PL+Un7/nnZyU0rFH00CgDGG4tLymOy7pKwCgNz9R2Oyf6WUqg5NAsAHi7bT4+Hv2L7vSNT3bbBmbhOJ+q5joqi4lH2HS5wOQylVQzQJAN+u2gVATsHhqO/bNXtnSoJkgdOf+p4Bf53maAyb9hTRcfwkZm7Y42gcSiUDTQIcv1qPhQo7CyRGCoBDx8qcDoElP+8HrAp1pVRsaRLwIDE8VUuC3AlEYlvBkZjUqZjY5WallE2TALE52Rhj+HbVLioqor9vf46VlfP81A0xq+T2Z/gzM7n530uitr9YJmOllDdNAhxPAtG8WP9q+U5ufX8pb87d4nff09fmkZ1XFLVj/mfBNl76fhOvzNoctX2GavbGfPdyaXkFz0/dwOE4KFpSSgWmScBDNK8/99qdw3YesJqG7vDRRPT37y1mxAs/RO2YrjuA0vLwbz8Wbd3HroOhNWOdm73XXV5vfNxGfbE0l5e+38Tz0zaGHYdSqmYl1Mxi0bJgSwFdshqS1agOEFnF8LS1eeQVFnPNaR0Cbrdw676g+9p18CitM+uFHUNlJoJK6Ce+Xkvzhhk8M2UD9dJTQ3rPNW8tBCBnwmifrx+zk1B1i6W0SkCp2EvKO4ErJy7g8ld/rNY+bnxvMQ99udrv63Oy9wJQVBy8SOQPEZSnHykp458zN1HmcdVfYZ81X5m1mdOfmhHSft6et5VnpmwA4GgEJ+1Ii9LyCovDPlY0rdpxMKpFcUolqqRMAgDbIuwYtnZnIec+NyvodnM37Q15n8GaZa7OPejVgaukrIJej0zhmSkbGPHCD/x34TbAu1x+18HYnGQ3VjpxRtIEdtaGPQx5cgbT1+ZFMbLw/OLluVEtilMqUSVtEvDkLtb2cyYrKatwX3E/P20jW/KDdypLCeOsGKzY4+L/N5dLXp7rfj7+s5Xu5a17D/PAF6sA2HUgdkNT/Lh5L1PW7Oa+T1b4fN1XE1h/n2vljoMArNhxwPcG2jhIqRrjWBIQkXYiMlNE1orIGhG5K9bH3LTnEC9Oz/Yfk5+zz0kPTWbMP+eFdaxoN3P0rFj+fFmuz212hnD1v2jrPncP6XBc/cZCn81Aq1NuP39zAbsDxKz9BJSKPSfvBMqAPxpjegGnAbeJSK9YHvDKiQt4YfrxFivv/phDcWl5SCeyNTsL/b5mjOGr5bkcKzteph5WGXkNnuzGvj6fW99fGvkOKn0wd52Ar02D7Grxz/u56MXqF8mUllewrSD64z4plQwcax1kjNkF7LKXi0RkHdAGWBurY5aUeVd8Pvq/Ne4mnBB5P4HZG/O568PljOl3ImMHtePMri04VlaDvcQcFOoAeatzDzJ1zW7uHdnda/3+I6UhHWfljgPM31zAzWd18Vrf9YFvKbNrxJc+PIJmDTJCjFwpBXFSJyAiHYH+wEIfr90kIotFZHF+fn6V94bD1wX3waOlfq/E/2/SWro/NDnofgvtFkBfLd/JuDcXsj/MUTidKvUorwj/yP7O9aXlVff1vl1hDVZF7Evfb4r4GJe8PI+nJq+vsr7M4zMUFYeWUJRSxzmeBESkIfAZcLcxpkqZizFmojFmkDFmUFZWVmxj8Viek53PG3O2hnRFf+cHy7yeh9J6qKaHdqhs1Y6DdHngW7+vf7NyJyt9VNwu3358XUWF8Sq3/8f0jRhjfJ7EPbcLdsNVYlfC++qIppSKLkeTgIikYyWA940xn8f8gD7OKb7OM2t3FnLt24tC2omvq+lQijimrNnt9byiwvDlslyv/RUcOkbH8ZOC7sufsvIKr3oKTz/lBO7Edvt/l3HJy1ZluL/K2+/W7K6UBLLZtOdQwP36KrsvrzBen9tVef+NjwrsZ6focBRKRZOTrYMEeAtYZ4x53qk4oGqP4VEvzSHUkpJIilTg+HDJYF3xfrpkB3d/tJx35m11r1+27UBE+3YZ+/p8uj/0XbX2sTn/kN9x/VflHqyyrtx4f5uVr+ZfnV21SKj/E1MZ8uTxzm37j1jFaSU+7sJenrmJF2f4b+GllAqPk3cCZwK/Ac4VkeX2Y5QTgbh69R44GkqZcqXWMWGU6P/pkxUcOlZGUXEp783/2eu1ArseIT+KE9IvDZBEQq24Pu+52VXuWlxKyyqqfP5ISnAKi8vcYy2Bd5+D+ZsL+LFSxzuni9KUqk0cSwLGmLnGGDHG9DHG9LMf/gupo3FMP+vX77Z6wb4a49E3P12yg3d/zOGHjd4ntZyCI2zbZ3VAW7Xj+NV15RY3i0IYh8iXjuMncbTE+8T5xpwtIb9/1gbfFfJlleoEwEoCnmF/s3IXuR4tsD5YtD3o8Tzff9UbC7j6zSrtBZRSUeJ4xXBNCjY8g5PzvqzOterEf9xcwMQfNtNx/KQqLW7Gvj6fwghbwBQcPuZVvBKNq+l//ZhTJbFWvjO444NlXP5KaOM0dRw/ifyiY0HvUoyBPYXF2hpIqShIqiTgy0eLg1+ZeiqpNExzuMUfFRWG1TurlqV7njxdFaOPf72mynZ9Hpvqc7+hVCD/bVL0u2BUrgj29X3sDmOwuM35gSuWXQY/OYNzn5vtte6Jr9dWGdtIKRVY0icBT5WLTHz5YaN30Ui4V9TPTdvos9jJ89iH7eVoDwLnWQ9xJITPGopLKw2nkV8UvToNf1wJs/KxZqzfw8hqDAr31fJcOo6fpHcYKqloEvDgqhsIx7golVdH66TsT001ub/+Xz/5nEAnVNGI03OY6py9h/l86Y6Q3vfKTCs558ZwID6l4k1STioTyN4wW+cEGlMoHAdDapkUOX8tfGIh2Axlh0qq184/WKIoPFpKq8Z1ARj90hwOl5Rz2YC21TqmUrWV3glUcv07Pzly3FjfCfxt0rqY7j8c+w75H1Yjklnequ7juMMx/l6VSnSaBCrRooDq+2r5zoCvB6rruPqN4MVroaaJb1YGjsPv/j0O8NXyXA6GOMidUolIk4CqceHMuuZLsOIg1+sb80JraeRSuYnwF8t2cNeHy7nzw2W+36BULaBJoBKd1CrxRaNICeCej6xZ1AJNfKNUokuaJLBgS4HTIagoCbnITkchVSqopGkddOXEBSFtF6jX8H89xsdXzqncV6OyQOf+jXlFpAaZALry+53sSa5UrCVNEgjV3gAtV1wTuqv45jqJ+8oF1elMplRtlDTFQUopparSJKBqnepWDGvxj0omWhykap1FW/exac+hkOuFp63No02Teu7nWp+skokmAVXrPP511dFSb3t/KZ1aNPC5/Y3vLQagxwmNqn3sA0dKKCmvoLzC0DqzXvA3KOUwTQIqKUzyMV+xP2/N3cpzY/tGdJx+T0xzLy984Dz3GEZKxSutE1Cqks8CjDq6ac8hvg0xoRQEaGmmVLzQOwGlwnD+89ZENjkTRjsciVLRoXcCKql5TiAjfpoF+Vq/Jf8Q5z43i4Iwhx5XKt5oElBJ7ab3lriXt3hMbbk6t+oUoJ4m/rCFLfmHmbo2L2axKVUTNAmopDbfY0wpzwnuXS2GPBkfbUe1OalKdJoElAqR58Q8rhKiaI1YqpRTNAkoFYSrRuCtuVt9rFUqsWkSUCqI4jL/U1QaA7M35rM5P7wJbJSKF9pEVCkfPK/zt+Qfrvq6uzgIrn17EaDNRlVicvROQEQuFJENIrJJRMY7GYtSnoKV9GthkKotHEsCIpIK/BO4COgFXCUivZyKRylPuypNKfni9GzfG2rzIJXgnLwTGAxsMsZsMcaUAB8CYxyMRym/Xpi+0eu5qziosLgspPd/uGgbW/dWLVZSymlOJoE2wHaP5zvsdV5E5CYRWSwii/PzA08rqFRNEbtA6JkpG/xu41mhPP7zVVzy8tyYx6VUuOK+dZAxZqIxZpAxZlBWVpbT4SgFQE5B8Kv6O/67zOt5UXEZT01ex77DOrCcih9OJoFcoJ3H87b2OqXi3pzsvUG3yT1wtMq612dv4fGv18QiJKUi4mQS+AnoJiKdRCQDuBL4n4PxKFUtB4+WBt8IKC2vYMf+Ixw+5r8+oay8gldnbabEYygLpWLBsSRgjCkDbgemAOuAj40xeomkEtYfP14R0naCMPTvM/n1xPnudYePlfHi9GzKyq2T/rNTN/L379bT9/GpMYlVKRdHO4sZY74FvnUyBqWiZfq68EYUXZ1b6F5+duoG3pmXQ9um9bh8YFt3D+Sjpf57KysVDXFfMRwNP24OXn6rVI3x0dPsaIl1si8p1+IfVbOSIgks2LLP6RCUCon2PVM1LSmSgFauqXjia8gJP5OaKRVzSZEEXpu92ekQlHL7ZqX3RPUzN+xhzU6rfqCwuJQfN2nxpao5OoqoUjG0blchPVs39vv6wi0FXP/OT+7nEyavB+CMLs1jHptSkCR3Ako55aIX53DwiP/+AwV+eg97TnVZVBxa/wOlIqFJQKkYe3rKer+v+asIXvLzfvdy78emsmzbft8bKlVNmgSUirG9h45Vex+rcg9GIRKlqtIkoJSDwp2W0hjDxB82s6ewOPjGSoVAk4BSDioMcbwhVwvS7D2HePLb9dz236WxC0olFU0CSsXYlDX+h5MIt29YWbn1jqIQJ7NRKhhNAko5KOQewnZvMvcE99qzWEWJJgGlHPT2vK0hbffwl6sB7Vmsok+TgFIJyIRdkKSUb5oElEoQxaXl7rmNlYoWTQJKJYipa/NYt8saY2hj3iH38NOVfbFsBx3HT4pK/wRV++nYQUoliDs/8J64vucj33HXed24Z8RJbN93hGFPz+Sqwe35YNE2AHL2HqZFwzpOhKoSiN4JKJXAXpyRzU85+xj29EwAdwJQKlSaBJRKcFvC7HWslCdNAkrVUtqcVIVCk4BSCe7Pn61yOgSVwDQJKKVUEkuKJNC2ab2Arzeqo42kVO10tKSch79crRPTKL+SIgns2H804OsXnnJCDUWiVM16d34O/17ws86zrfxKiiTgzyvjBgDQXNtSq1pJKK+whpeo0FEmlB+OJAEReUZE1ovIShH5QkSaOBHHRaecwITLenPPiG5R3/c1p7WP+j6VCkdJWQWvzbLuAMIZdTRPJ6xJKk7dCUwDTjHG9AE2An9xIggR4crB7amTlhr1fd80rIt7eUinZlHfvz9nd8+qsWOp+Hb3R8soOhbevANLt+1nyJMz+HTJjhhFpeJNSElARO4QkabROqgxZqoxxvXrXAC0jda+o+22c7oE38iH9s3ru5dH9GoVrXCCevjiXjV2LBXf8gqPjx0U6qijG3YXAfDT1n0xiUnFn1DvBFoBP4nIxyJyoUhUu6HcAEyO4v6iSkdtVLXZ0ZJy5mTnu5/rrz35hJQEjDEPAd2At4DrgGwReVJE/F4mi8h0EVnt4zHGY5sHgTLg/QD7uUlEFovI4vz8fH+bxZUzuzb3eh6oPPYPZ0V2p6FUdYz/bCUdx0/iwS9W8Zu3FlWZ8F7nK0geIdcJGGMMsNt+lAFNgU9F5Gk/259vjDnFx+MrABG5DrgYGGfv299xJxpjBhljBmVlJW5599e3D/W5/t4RJzFv/Lk1HI1KNl8t2+n1/MOftgOwyT75H7LnLNbpK5NPSL2kROQu4LfAXuBN4D5jTKmIpADZwP3hHFRELrTfc5Yx5kh4IdesaFwRGQy922ZWWT/n/nPISEuhTZPAndlCded50W/lpGqH3driR/kR6p1AM+AyY8wFxphPjDGlAMaYCqyr+XC9DDQCponIchF5LYJ9hKxOWuSNoDJS/bccGjfkeDPQtBRhVG+r01nnFg29tvN3VdWuWX3fL0To3hEnRXV/Knms313Ixrwip8NQDgh1vIQXAUTEs61jkTGm1BizLtyDGmO6hvue6vjVoLb8Z0Fk46zfNLwzL0zf6H5eNz2F4tIK+rVrwmUD2vD+Qmu/m54cBcAPG/MZ0tl3k9AuWQ3YnH84ojhClZ4SWcL7+vah/OLluVGORsU71wVK5UHotDQoeYR6xlgK5GO16c+2l3NEZKmIDIxVcE5KTRGaNcigXkYqf/S4wr7ujE4ApKf6bkcx/KSsKv0OXP+h0lNj3y3Ds2lqqHImjPZZXKVqpwVbCpwOQcWRUM9K04BRxpgWxpjmwEXAN8CtwCuxCi5W6qSl0DqzbsBtpt97FksfHgF4j8vuWUcQrPLst6d38Hr+7K/6hhdohM7RDmMqgCsnLgi6jTFw3ycr6Dh+Ug1EpJwUanHQacaYG11PjDFTReRZY8zNIpJwA++sfeLCoNukpRw/80faUqJlI+urqZ9h3Rmc0iaTnAmjOXiklGNlvicJd/n37wbzm7cWRXZgpXy475MVPFPpQiRQw4dPtNdwUgj1TmCXiPxZRDrYj/uBPBFJBSpiGF9MpKYIqSnexTldshoA8NDonozs1cpvpW0Xu9L36iHtqZ8ROIfeOLwzD47qydWDvccRyqyfTsvGge9EhnUL72q+b7sm7uVo9eV7aHTPkLY7t0fLqBxPxdYnS3bolb2qItQkcDXW0A5fAl8A7ex1qcDYmEQWRZ5X8q0a+75x+fK2M/nhvnP4/bDOTPztIK/XPM+pzRpkkDNhNL/s39adSE5q5d0ayKVOWio3Du9MWhTqAr65w3c/A7CKf/7zu8FV1vs6if+yf5uQjte9VSN+P6wzORNGB9yuV+vG3KVNU2sd7SyWPIKeneyr/ReNMXcYY/obYwbYy/nGmBJjzKYaiLNaXCfxeumpzL7vHJ/bNKqbHlKlqvFaju1/FM8T8Clt/Ffcds5qSKO66e7nrpzVsXmDKts+PzZ4vcTXtw/lo5tPCynG83u1om+7JkGThUowmgOSRtAkYIwpBzqISEYNxBMTKXYWqJueQt306I8Y6pQ7zo1NS9vebTNpUj+0P7eONZPYVucW+lyvOSB5hFoxvAWYJyL/A9wN3Y0xz8ckqihzJYFIJ9YY2i2LZ6dafQUCjHARFVPuHk6DOqlez0vLfVe7DOxgDezaq3Vjn697Rjp2UFvaN6sfsL7gmzuG8nNBdDpwZ6SlUFIWuLqoT9tMVu44GJXjKaUiE2oS2Gw/UrB6+iYU13mvIsIs0M+j0tXn/qNwPXzr2V1Yt6uQ7id4f72Vn3s6u3tLpt0znK4tveskfJ3nn76iajHQU5f1ZlCH4yOEn9Im02ex029P78AZXVowY10elw1oy1VvHG9i6HmsK09tx7S1ebx/4xCa1c9g8JMz/MYO8L/bh2pFZQI6eLSUz5bsYE52Pu9cX7UuSiWWkJKAMeZxABGpH+9j/fiS6r4TiPwq/vyerZi+Ls9rXTRvCu6/sEfQbZ77VV9OatWIS1+Z5542sFuryHPyVYNDm/3siTGnANZczMu27fe73YTL+zDh8sD7euaKPtz36cqQY1TO+GJZrt/Xzn9+NvlFx/y+rhJLqAPInY41jHRDoL2I9AVuNsbcGsvgoqVxPavSdHSf1tXel5NlpZcPtObeWfHoSHcSCCRQ0VV1v4s2TerROrMuVw/xn0iGdWvBnOy9XuuiOxWFcoImgNol1LaL/wAuAAoAjDErgOExiinqXJ21HhgVWrt3X+Lp3NWwThqZ9dIDbHE82Awfg+cteeh8Xhjbr1oxtGhUh09vOYOWjfz3d/CVg6rzNY4LkHA83XVeN1Y8MlKbrkbZwaOlFBaXOh2GirJQ6wQwxmyvdBUXuMtrHEpJqf6Z3OeJLY4SBMD4i3pQeLSUM7u2YO7957D3UInX680bRt7J++QTMzm/Z0v+OLJ7RO8PVMcRTKjf8z32WE/3jDiJF2dkAzCyVyumrs0L9DYVRN/HpzodgoqBUO8EtovIGYARkXQR+RMQ9uihicxVOew59n+8TrzRtWVDPv7D6TSok0bLxnXpdaLv1kORyEhL4c1rT6WnnxZJnsae2s7r+a1nd+GUNpnU82im292u0/A8wU+5ezhPXdYbgKxGCTcqiVIJJdQk8AfgNqANkAv0s58nlOpcsN9yVhem3D1cR9sMwyV9T/TqRNasgdX34Lu7h/Hy1f0B+OAmq1Na5xbHO7Z1P6ERJ9gD/Hk2f63cCuv+C4Pfjbx4Zb8aG7ivtvr9uz85HYKKoVBbB+0FxsU4lpiJxhV7SopUKcronNWA7q0a8dglJ1f/ALVYaopQXmFonWndRXVo3oAOdm/mZg0yePeGwfRq3ZhH/7eauvYw3Kd3bs7IXq34y6ienPPsLKDqXcGtZ3fl6e82BDz2mH7WMBlT1uyO5kdKKtPX7dGmvLVYqK2DsoAbgY6e7zHG3BCbsGIj2i1T6qanMuWehKkfd8zZJ2UxY/0en5XUAGedZA2W98q441NT1E1PrTKG0y1nd6Ft03qM6t06bovilEo0oVYMfwXMAaaTgBXCylmu83V1U3B6agqXDWgbWQyaNJTyKdQkUN8Y8+eYRhJDOiKis1z9FZxsRXVa52ZVOvup6jt4pJTM+oGaK6t4F2rF8DciMiqmkdSAOGvJmTSuGGi1EgqlRVGs/G5oJ7+vDe7ke05oFdzk1bu8nq/YfiDohEkqvoSaBO4CvhaRoyJSKCJFIuJ7+ME4pEUBzhrdpzU5E0Zzokfz2mhZ+8QF9G/fJOh2IuKzd3Pfdk144zeDmPmns6MeW7LZvu8IY/45j0e/WuN0KCoMoSaBTOA64CljTGPgZGBErIKKlXjr1KWqr35Gmle/g0Ce/GVvzu/pPQvaV7edSWb9dDq1qDr3QmXr/xp8WtJkU1hcykJ74voDR6zexKt36siwiSTUJPBP4DTgKvt5EfByTCKKAb0RSGy/Oa1DwNdvO6croXYGD3RXeP2ZHQO+t256Ko3qhNzJPik8+e16fj1xAUU6nETCCjUJDDHG3AYUAxhj9gMJN8lMNIZ8VjXvr5eeEnDmsjO7tmDLU9Wf2ezRX5zMikdHBtzm6tOqFilFYTSShLcx71CVdQeOlPicC2P59gO8Nz+nBqJSoQg1CZTa00wacPcbSJgJ5rVOQIUqs146ORNGM7JXK98b+Pgt6dAWkJ1X5F52/X/r98Q07v14RZVtL/3nPB7ReoO4EWoSeAlrgvmWIvJ/wFzgyZhFFSNaJ6BC9eKV/Zl81zDuPt8aibRNlCq1A7VSSmSzN+b7bIr99YqdDkSjwhFSEjDGvA/cDzwF7AIuNcZ8Ut2Di8gfRcSISIvq7kupUIQ6d3K9jFR6tm7MsG7WT7NVY+tq3/M05+rpfN8FwScEcrlqcLvgGyWgyat38/6Cbe7nS37e516+/p1FLPnZ/2REylnhDCW9HlgfrQOLSDtgJLAt2LbVpZ3FksNb1w6iR5C+CI+POZk+bTN59H/VL47ISEtx11X86ZOqxR6+1OaiyY8WbwesO+4//Gepe/3MDflk7znE3D+f61RoKgAnmzq8gHV38ZWDMaha5LyefsrxPTSsk8a1Z3SkSf10WjX2PyGOS6CTtudr/ds3Ydm2A1W2OalVQ69K0zppoTVnTWSrcwu1niSBhFonEFUiMgbItWcoi7nafPWlIjOmXxtO69w85O2DDT44oH1TAFpneieWK+wpQTNSU3h13ADaN6/PeT1aVnm/S4uGCdfoTiW4mCUBEZkuIqt9PMYADwCPhLifm0RksYgszs/Pr2ZM1Xq7SkL1M6ybZdfJ3d+8za7VvxvaiSEew1C4miU3rpfGRb2teZ3fuu5U99wKvxroPSDeiF4nRC/4OKL/9+JXzJKAMeZ8Y8wplR/AFqATsEJEcoC2wFIR8fnrN8ZMNMYMMsYMysrKilW4SvnU68TGvHRVfyZc3ifgdv7qnVwnv8rFQJf0PRGgSh3G0K7aRkLVrBqvEzDGrALc98N2IhhkT1wTU9pZTEXCdcIOxHUnICJe6aB5wwz+NPIkRvfxvY/a+ovMLzrm9dzf/71Dx8pYvu0AQ7tp8nOKI3UCNc3fLbxS4fL8Kd16Thf3ct921rSjrjmTPd1+breQxiaqzbbtO+Jz/b0fLeeatxay88DRGo5IuTg+EIoxpmNNHUvLJVU0uSqDAX7Zvy0D2jelQ/MGvPR9tnu93n36d/BIKZv2WC2njpbq8NNOSYo7AaWiJdA9pWve5FNOzAy+nziYaKemXfHqj17PCz0GndObdeckRRLQH5iqSX8Z1SNo8U+kU27ecnaX4BvFqcWVeg0bQ+2tFEkgjhcH1ST9vanqcl1QjOnnv7I4PTWFPm0z2br3cND9Bet/AHD3+d3o2rIhDeukcXb3lrw6a3PI8ca7Y6WucSj1Ss0pSZEE9Oeloq13m8BFPsHuPv293rqJ717MF/tpXZTI1uw8SK5WCDsuKYqDXEK56lIqkF8Nsjp3jaxmp64uWVZxUdum3qOTelY213bLdxxwL2uRrXOSIgnoD0xFS8/WjcmZMJr2zetXaz+/Pb0jH998Ouf1bMXku4ZFtI8Jl/Xm5BMDD5gXb979McfpEFQlSZEEXPQ+QNU0fzefKSnCYHt4iZ6Veg1fGqC+wdOVg9tzwcmJNcyEv9Fb9SbdOUmRBHQoaVXTqvOL+8eV/Vn+yAjG2kVPJ4fQ5DTR6d26c5KiYthFrzZUomhSP4Onr+jLHed2o12zqkVPVw2uOtdxotGOdPEhqZKAUonGVwJwTWRTm+iNgHOSozhIf2FKxb3dB4vJKyx2Ooykk1R3AtpEVNU2vn7RIol34WMMnPbUDKB23unEs+S4E3A6AJV0Hh7dk8v6t3Gk9c6yh0fU+DEjoddk8SEpksBLM7KDb6RUFLVsXJfnf92Puuk1M6fw7ed0dS83qa9TVKrQJUUSUEop5ZsmAaVqgUTsC3PgSIl7uaSsIsCWKpY0CSiVRK6oNLG9kz5YtN29/OzUDQ5Gktw0CSgVBxrXTaqGelVo01DnJPcvT6k4Mf8v51FWEb0inQHtm7B024Eq6+ukxed13/rdRU6HkLTi8xehVJJpUCeNzHrpUdvfZ7ec4V4+tePx4an/fFGPqB1D1Q56J6BUAvPX1l5E6HFCIw4dK+Ojm07n9R+2sDhnH43rptO2aT127I/fyVzmbdpL7oGjtG1ajzO6tHA6nFpPk4BStdR3dw93L1tzE1vzE3908+mcOeF7h6IKbtybC93L2ns49rQ4SKkE1rKRNR1lq8a+p6X05cTM0LdVtZ8mAaUS2BUD2/LKuAFcM6SD06HExLaCI+wpKmZu9l6nQ6m1tDhIqQSWkiKM6t3a6TBiZvgzM2nTpB65B45q0VCMaBJQqpZY9OB5lJYnXs/hYHIPxG8ldm3gWHGQiNwhIutFZI2IPO1UHErVFi0b1aVNk3pOhxEza3YedDqEWsmRJCAi5wBjgL7GmJOBZ52IQymVOB74YrXTIdRKTt0J3AJMMMYcAzDG7HEoDqWSjojw98t78+VtZ3qtn3TnUK4eEr9zF+v0A7HhVBI4CRgmIgtFZLaInOpQHEolpV+f2p5+7Zow5/5z3OtOPjGT+y/o7mBUgS3ffoDsvCIqoji8hophEhCR6SKy2sdjDFaFdDPgNOA+4GPxM/ejiNwkIotFZHF+fn6swlUqKbVrVp9v7hjKW9cO8lrfuG4aORNG88jFvRyKzLcRL/zAyzM3OR1GrRKzJGCMOd8Yc4qPx1fADuBzY1kEVAA++4cbYyYaYwYZYwZlZWXFKlylktYpbTI5r2crAKRSocsNQzs5EVJAy7btdzqEWsWp4qAvgXMAROQkIAPQ3iBKxYl4LnDxU2igIuRUEngb6Cwiq4EPgWuNMTH73V02oE2sdq1U7aLn16TjSGcxY0wJcE1NHS+zXjqNknzSDqVqiw0690BU6dhBSim3RnXSGNC+CS9e2c/pUPxy9SD+17yt/Fxw2OFoEp8mAaWUW0qK8PmtZ3Juj1budfFYBH/gSAmPfb2WKycucDqUhKdJQCkVUKfmDaqsG9ihqY8ta85fPl8FwK6DOjdxdSVFEigureDwsTKnw1AqMfm4E3D65mDy6t0OR1B7JEUS+GDRNrSToVKRubjPiQA8fXkfhyPxraSswr1cXFrOyh0HnAsmASVFElBKRe7u87qx8rGRjD21ndOh+PTm3C3u5Ye+XM0lL8/T4afDoElAKRVQSorQuG6602H4dbSk3L28YvsBACav2qVFwCHSxvNKqYSWIsKSn/excOs+Kuw+p3+btI6l2/bzyriBDkcX/zQJKKWion2z+mzbd6TGj/vijGxenJENQOes4y2ZtuQf70NQWl7B8u0HOLVjsxqPL95pcZBSqlpc/QhMHIw45Hni9xxj6NmpG/jVa/O10tgHTQJKqbB5diBzurmoPykega3fZQ01UXCoxKFo4ldSJIGhXX2OUq2UCtNd53Xjpav688ConnRoXp/mDTJ4+eoBANx2dleHo/O2Zmehe5wh1z1KPPZ+dprEcPDOqBs0aJBZvHhx2O8rKauguKw8rls4KFUbXPLyXFbuiL8J4Vs1rkNe4TF6tm7MO9edygmZdQNuf+E/fqBt03q8eW3tmPRQRJYYYwb5ei0p7gQy0lI0ASiVxPIKjwGwblch9368POj263cXMX1dckx9nhRJQCnlnHgrgim3hw847ckZvD57s8PROE+TgFIqanyVLo/u3brmAwkgRYSFWwrYXVjMU5PXOx2O47SfgFIqpuqmpzodgpf5WwqYP7HA6TDiht4JKKVUEtMkoJSKmocv7uX3tR4nNKrBSFSoNAkopaJmcKfjwzK4Omu56oWvPaNjjccTikPHyvjzpyspKi51OhRHaBJQSsXE0G5ZADSoY1U91kmLz9PNXz5fxUeLt3P9Oz85HYoj4vOvopRKWPXsiuDXrxnI9388i/sv7M59F3Tnkr4nOhyZb7n7rUHvFv+8n8e/XhN0+3/Pz+Gv36yNdVg1RpOAUiqqpt4znHeuO5V6Gal0zmpI/Yw0bjunK2mp8Xm6WbrtgHv5nXk5Qbd/+Ks1vDV3a+wCqmHaRFQpFVXtmtWnXbP6Pl9rVCeNojif7KXj+ElOh1Cj4jM1K6VqpUcvOZlGdRLn2rPj+Ems313oNY+x52u1gSYBpVSNuWJgW1Y9fkGV9S0aZjgQTWgu/MccBv5tGh/9tK3Kib+iwv8AnJ8t2UF+0TG/rxcWlwZ8vaY4kgREpJ+ILBCR5SKyWEQGOxGHUsp5H998Ot/eOYym9eN3kMei4jL+/NmqKuvf+TGHB79YVWU+4z2FxfzxkxX8/j3/ox73eWwqp/7f9KjHGi6n7sueBh43xkwWkVH287MdikUp5SBX34JGddPZfySx2uo/9e06yioMJzapx23nWPMpdBw/iTO7Ngcgv7DYyfBC4lRxkAEa28uZwE6H4lBKxYmKBJrbxKXMT3HQvE3W2ESJ8ImcSgJ3A8+IyHbgWeAvDsWhlIoT/nLAkofOr9lAIlBUXMbhY2VVeh3vOljMNW8u5POlO9hTVMzbc7eyOje+Jt2JWXGQiEwHTvDx0oPAecA9xpjPRGQs8Bbg8y8tIjcBNwG0b98+RtEqpZyw8W8XuZebNcgg98BRnhhzMnsPlTBuSHvqpqWSWT+dC05uxZQ1eQ5GGthrszfzmp+5CeZu2svcTXvp2zaTFfasa29fd3ySr1v+s4RxQzowtJsz0+A6Mr2kiBwEmhhjjIgIcNAY0zjY+yKdXlIpFV9crWxyJox2r8srLGb6ujzGDelQZfuJP2zmyW9r/9j/H998On3aZlI3PZU352zh3B4taZ1ZD4B6GZEPyR1oekmnKoZ3AmcBs4BzgWyH4lBKOSSj0lhCrRrX9ZkAAG4c1pnzerai4FAJY1+f7/XasG4tmJO9N2Zx1qSxr89nYIemPPnL3vxt0jpenbWZgsMlpAhseWp08B1EwKkkcCPwooikAcXYxT1KqeTw8MW9GB5G8YeI0CWrIV2yvNePHdSWe0d054EvVvH9+toxJ/CSn/dzwT9+AKDgcAkAAbojVJsjFcPGmLnGmIHGmL7GmCHGmCVOxKGUcsbvhnaiW6vI5hf47JbTAejZujFPX9GXEzLr8vpvBvLkL3u7t/n1oHZRiTOezNwQmySnPYaVUgllYIdmLH7ofD6/5Qz3uvTUFK4e0p7+7Ztw+zld+fsVfVj+yAh6nNCI6+J0HoNwzVgXm4rxxBnEQymlbC0a1vG5/otbz3QvN6mfwXd3DwfgXz/mVPuYN5zZibfnRW/00OWPjODzpbl8uTyXnic0pqzC8NnSHX63z9l7JGrH9qRJQClV6215chRLt+2nfkYaBYePcUaXFtz87yVs2XuIe0ecxGmdm9OkXjpT1uTx7o85XNT7BB7/ei1pKcI3dw7lp5z9XDOkPR2a1+ezpTsoOFRC7oGjfo9XLz2Vo6Xlfl/f+tQoRIQbhnbihqGd3OsDJYG/XnpKZB8+CEeaiEZKm4gqpVT4AjUR1ToBpZRKYpoElFIqiWkSUEqpJKZJQCmlkpgmAaWUSmKaBJRSKolpElBKqSSmSUAppZJYQnUWE5F84OcI394CiMfxZjWu8Ghc4dG4whOvcUH1YutgjMny9UJCJYHqEJHF/nrMOUnjCo/GFR6NKzzxGhfELjYtDlJKqSSmSUAppZJYMiWBiU4H4IfGFR6NKzwaV3jiNS6IUWxJUyeglFKqqmS6E1BKKVWJJgGllEpiSZEERORCEdkgIptEZHyMj9VORGaKyFoRWSMid9nrHxORXBFZbj9GebznL3ZsG0TkgljGLSI5IrLKjmGxva6ZiEwTkWz736b2ehGRl+zjrxSRAR77udbePltErq1GPN09vpPlIlIoInc79X2JyNsiskdEVnusi9r3IyID7e9/k/1eqUZcz4jIevvYX4hIE3t9RxE56vHdvRbs+P4+Y4RxRe1vJyKdRGShvf4jEcmoRlwfecSUIyLLHfi+/J0fnPuNGWNq9QNIBTYDnYEMYAXQK4bHaw0MsJcbARuBXsBjwJ98bN/LjqkO0MmONTVWcQM5QItK654GxtvL44G/28ujgMmAAKcBC+31zYAt9r9N7eWmUfpb7QY6OPV9AcOBAcDqWHw/wCJ7W7Hfe1E14hoJpNnLf/eIq6PndpX24/P4/j5jhHFF7W8HfAxcaS+/BtwSaVyVXn8OeMSB78vf+cGx31gy3AkMBjYZY7YYY0qAD4ExsTqYMWaXMWapvVwErAPaBHjLGOBDY8wxY8xWYJMdc03GPQZ4115+F7jUY/17xrIAaCIirYELgGnGmH3GmP3ANODCKMRxHrDZGBOoV3hMvy9jzA/APh/HrPb3Y7/W2BizwFj/W9/z2FfYcRljphpjyuynC4C2gfYR5Pj+PmPYcQUQ1t/OvoI9F/g0mnHZ+x0LfBBoHzH6vvydHxz7jSVDEmgDbPd4voPAJ+WoEZGOQH9gob3qdvuW7m2P20d/8cUqbgNMFZElInKTva6VMWaXvbwbaOVQbFfi/R8zHr4viN7308ZejkWMN2Bd9bl0EpFlIjJbRIZ5xOvv+P4+Y6Si8bdrDhzwSHTR+r6GAXnGmGyPdTX+fVU6Pzj2G0uGJOAIEWkIfAbcbYwpBF4FugD9gF1Yt6NOGGqMGQBcBNwmIsM9X7SvHmq83bBd1nsJ8Im9Kl6+Ly9OfT+BiMiDQBnwvr1qF9DeGNMfuBf4r4g0DnV/UfiMcfm383AV3hcbNf59+Tg/VGt/1ZEMSSAXaOfxvK29LmZEJB3rD/y+MeZzAGNMnjGm3BhTAbyBdQscKL6YxG2MybX/3QN8YceRZ99Gum6B9zgQ20XAUmNMnh1fXHxftmh9P7l4F9lUO0YRuQ64GBhnnzywi1sK7OUlWOXtJwU5vr/PGLYo/u0KsIo/0nzEGxF7X5cBH3nEW6Pfl6/zQ4D9xf43FkplRiI/gDSsSpNOHK90OjmGxxOscrh/VFrf2mP5HqyyUYCT8a4s24JVURb1uIEGQCOP5R+xyvKfwbtS6ml7eTTelVKLzPFKqa1YFVJN7eVm1YztQ+D6ePi+qFRRGM3vh6qVdqOqEdeFwFogq9J2WUCqvdwZ6yQQ8Pj+PmOEcUXtb4d1Z+hZMXxrpHF5fGeznfq+8H9+cOw3FpMTYbw9sGrYN2Jl+AdjfKyhWLdyK4Hl9mMU8G9glb3+f5X+ozxox7YBj5r8aMdt/8BX2I81rn1ilb3OALKB6R4/JgH+aR9/FTDIY183YFXsbcLj5B1hXA2wrvoyPdY58n1hFRPsAkqxylN/F83vBxgErLbf8zJ2r/0I49qEVS7s+p29Zm97uf33XQ4sBX4R7Pj+PmOEcUXtb2f/ZhfZn/UToE6kcdnr/wX8odK2Nfl9+Ts/OPYb02EjlFIqiSVDnYBSSik/NAkopVQS0ySglFJJTJOAUkolMU0CSimVxDQJqLgnIkZE/uPxPE1E8kXkmwj310REbo3gfQ1F5HUR2WwPuzFLRIaEuY9ZIlKtycJF5FIR6VWdfSjloklAJYLDwCkiUs9+PoLq9RxtAoSdBIA3sQYl62aMGQhcD7QI9c0ikhrBMX25FGvkSaWqTZOAShTfYvWehEpjv9hjsX9pD1i2QET62OsfswcwmyUiW0TkTvstE4Au9tjxz9jb3iciP9n7eLzywUWkCzAEeMhYwyFgjNlqjJlkv/6lfXewxmNgPkTkkIg8JyIrgNMr7fMqe9z31SLyd18fWkQmiDX2/EoReVZEzsAaY+kZO/4u9uM7+/hzRKSH/d5/ichrIrJYRDaKyMXhfukqCVSnp6c+9FETD+AQ0AdrSOG6WL0szwa+sV//f8Cj9vK5wHJ7+TGsoTHqYF2xFwDpVB3mYCTWJN6CdWH0DTC8UgyXAF8EiNHVw7MeVm/N5vZzA4z12G4WVo/OE4FtWEMWpAHfA5dW2mdzrJ61rk6dTex//wVc4bHdDKy7E7AS1fce231nf6ZuWD1n6zr999RHfD1cAzMpFdeMMSvtoXevwror8DQUq+s/xpjvRaS5xyiQk4wxx4BjIrIH30P+jrQfy+znDbFOmj+EEeKdIvJLe7md/f4CoBxrsLDKTgVmGWPyAUTkfayJUL702OYgUAy8Zdd/VKkDsUejPAP4RI5PIFXHY5OPjXXnki0iW4AeWElUKQBNAiqh/A94FusuoHmI7znmsVyO79+8AE8ZY14PsJ81QF8RSTXGlHu9WeRs4HzgdGPMERGZhXXHAlBceftQGWPKRGQw1mQ7VwC3Y93peErBGnO/n7/dBHmukpzWCahE8jbwuDFmVaX1c4Bx4D4h7zWVxmivpAhraj+XKcAN9lU1ItJGRFp6vsEYsxlYDDwu4p5ntqOIjAYygf12AuiBNYJjMIuAs0SkhV1hfBUw23MDO55MY8y3WKNx9q0cv/05t4rIr+z3iIj09djNr0Qkxa7T6IxVvKSUm94JqIRhjNkBvOTjpceAt0VkJXAEuDbIfgpEZJ5Yk5BPNsbcJyI9gfn2+f0QcA1Vx4j/PdYEKZtE5CiwF7gPa0TIP4jIOqyT7IIQPssusSZUn4l1JzLJGPNVpc0aAV+JSF17m3vt9R8Cb9gV3VdgJcBXReQhrDqPD7FGigWr3mER0Bhr9MziYLGp5KKjiCpVS4nIv7Aqzz8Ntq1KXlocpJRSSUzvBJRSKonpnYBSSiUxTQJKKZXENAkopVQS0ySglFJJTJOAUkolsf8Pg4HQzCCWTc8AAAAASUVORK5CYII=\n", "text/plain": [ "
" ] }, "metadata": { "needs_background": "light" }, "output_type": "display_data" } ], "source": [ "%matplotlib inline\n", "import matplotlib.pyplot as plt\n", "plt.plot(range(len(energies)), energies)\n", "plt.xlabel('Monte Carlo step')\n", "plt.ylabel('energy')\n", "plt.show()" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "We can see that the energy is gradually getting lower and lower as the annealing progresses. It is useful when we want to know how to system looks like while Algorighm is running.\n", "\n", "## Get result -Result-\n", "\n", "We get spin sequence of calculation result using `result.get_solution`. In the case of classical ising model, we get it from `mysystem.spin` directly. However, `result.get_solution` is a useful method to get spin sequence from other systems." ] }, { "cell_type": "code", "execution_count": 13, "metadata": {}, "outputs": [ { "name": "stdout", "output_type": "stream", "text": [ "[1, -1, 1, 1, 1, -1, 1, 1, 1, 1, 1, 1, 1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, -1, 1, 1, 1, -1, 1, -1, 1, 1, 1, 1, 1, -1, 1, -1, 1, -1, -1, -1, 1, 1, 1, 1, -1, 1, 1, 1, 1, 1, -1, -1, -1, -1, 1, 1, -1, -1, 1, 1, -1, -1, 1, 1, 1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, -1, 1, -1, -1, 1, 1, 1, -1, 1, 1, -1, -1, -1, -1, -1, 1, -1, 1, -1, 1, 1]\n" ] } ], "source": [ "import cxxjij.result as R\n", "print(R.get_solution(mysystem))" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "This sequence is result of execution. It is expected to be in Hamiltonian groud state (or close to ground state)" ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "## C++ core interface\n", "\n", "We can do almost the same things with C++ core interface as above, with a few differences. We note two points the following. \n", "\n", "- need to input random number generator (C++11 random) to argument of seed.\n", "- in the Graph class, the access to $J_{ij}, h_i$ is slightly different.\n", "\n", "The content so far can be described in C++ core interface as follows." ] }, { "cell_type": "markdown", "metadata": {}, "source": [ "```cpp\n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "#include \n", "\n", "#include \n", "\n", "using namespace openjij;\n", "\n", "int main(void){\n", "\n", " //generate dense graph with size N=100\n", " constexpr std::size_t N = 100;\n", " auto dense = graph::Dense(N);\n", "\n", " //generate random engine\n", " auto rand_engine = std::mt19937(0x1234);\n", " //of course you can specify another random generator that is compatible with C++ random generator, say utility::Xorshift,\n", " //auto rand_engine = utility::Xorshift(0x1234);\n", " \n", " //Gaussian distribution\n", " auto gauss = std::normal_distribution<>{0, 1};\n", "\n", " //set interactions\n", " for(std::size_t i=0; i\n", " //std::vector> schedule_list;\n", " //utility::Schedule schedule;\n", " //schedule.updater_parameter = {0.01};\n", " //schedule.one_mc_step = 10; //number of Monte Carlo step per temperature\n", " //schedule_list.push_back(schedule);\n", " //\n", " //schedule.updater_parameter = {10};\n", " //schedule.one_mc_step = 80; //number of Monte Carlo step per temperature\n", " //schedule_list.push_back(schedule);\n", " //\n", " //schedule.updater_parameter = {0.1};\n", " //schedule.one_mc_step = 30; //number of Monte Carlo step per temperature\n", " //schedule_list.push_back(schedule); //schedule_list -> [(0.01, 10), (10, 80), (0.1, 30)]\n", "\n", "\n", " //do annealing (updater: SingleSpinFlip)\n", " algorithm::Algorithm::run(system, rand_engine, schedule_list);\n", "\n", " //show spins\n", " std::cout << \"The result spins are [\";\n", " for(auto&& elem : result::get_solution(system)){\n", " std::cout << elem << \" \";\n", " }\n", "\n", " std::cout << \"]\" << std::endl;\n", "}\n", "\n", "```" ] } ], "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 }