{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# A4 - 量子イジング模型"
   ]
  },
  {
   "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/A004-QuantumSystem.ipynb)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {
    "colab_type": "text",
    "id": "view-in-github"
   },
   "source": [
    "この章では量子 (主に横磁場)効果の入ったイジング模型をご紹介します。  \n",
    "まずはGraphを定義し、$J_{ij}, h_i$を決定しておきましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cxxjij.graph as G\n",
    "#問題サイズを100とします。\n",
    "N = 100\n",
    "\n",
    "graph = G.Dense(N)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "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",
    "        #Jijの値が大きくなりすぎてしまうので、全体の係数を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": [
    "グラフの設定方法は前章と同じです。"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## システムの設定 横磁場イジング模型\n",
    "\n",
    "今回はシステムに横磁場イジング模型\n",
    "\n",
    "$$\\begin{align}\n",
    "H &= s \\left(\\sum_{i<j}J_{ij}\\sigma_i^z \\sigma_j^z + \\sum_{i=1}^{N}h_i \\sigma_i^z \\right) - \\sum_{i=1}\\Gamma (1-s) \\sigma_i^x \\\\\n",
    "\\end{align}$$\n",
    "\n",
    "を用います。\n",
    "\n",
    "$\\Gamma$は固定されたまま、$\\beta$、$s$を変化させて量子モンテカルロ法を行います。\n",
    "デフォルトでは鈴木・トロッター分解による量子モンテカルロ法が実装されています。\n",
    "\n",
    "> 連続虚時間量子モンテカルロ法も用意してはいますが、現在試験的実装となっています。\n",
    "\n",
    "まずはシステムを生成してみます。`system.make_transverse_ising`で生成できます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cxxjij.system as S\n",
    "\n",
    "mysystem = S.make_transverse_ising(graph.gen_spin(), graph, 1.0, 4)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ここで、1つ目の引数にはスピン列を、2つ目にはグラフ、3つ目には$\\Gamma$の値、4つ目にはtrotterスライスの数を入力します。\n",
    "これで、全てのtrotterスライスが `graph.gen_spin()`で初期化されたシステムができます。\n",
    "\n",
    "`mysystem.trotter_spins`で全てのtrotterスピンを表示します。縦方向が空間方向、横方向がtrotter方向です。\n",
    "全てのtrotterスライスが同じスピンで初期化されていることがわかります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[[ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [ 1.  1.  1.  1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [-1. -1. -1. -1.]\n",
      " [ 1.  1.  1.  1.]]\n"
     ]
    }
   ],
   "source": [
    "print(mysystem.trotter_spins)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "> `graph.gen_spin()`の代わりに上の二重Listを入れて直接trotterスピンを初期化することができます。\n",
    "\n",
    "## アルゴリズムの実行 -Updater, Algorithm-\n",
    "\n",
    "### Updater\n",
    "\n",
    "量子モンテカルロ法に対しては、現状\n",
    "\n",
    "* SingleSpinFlip (メトロポリス・ヘイスティング法によるスピン1つずつのアップデート)\n",
    "\n",
    "が使用可能です。\n",
    "\n",
    "### Algorithm\n",
    "\n",
    "#### スケジュールリスト\n",
    "\n",
    "スケジュールリストは`(パラメータ, モンテカルロステップ数)`のリストで与えられ、横磁場イジングモデルに対しては(($\\beta$, $s$), モンテカルロステップ数)で与えます。例として以下のように設定してみましょう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "schedule_list = [((10, 0.1), 10),((12, 0.3), 80),((10, 0.8), 30)]"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "この場合、逆温度$\\beta=10, s=0.1$で10モンテカルロステップ、$\\beta=12, s=0.3$で80ステップ、$\\beta=0.1, s=0.8$で30ステップの計120モンテカルロステップを実行することを意味します。  \n",
    "アニーリングを実行するにあたっては、以下のように`utility`にある`make_transeverse_field_schedule_list`を使うとより便利です。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "[((beta: 10.000000, s: 0.000000) mcs: 20), ((beta: 10.000000, s: 0.111111) mcs: 20), ((beta: 10.000000, s: 0.222222) mcs: 20), ((beta: 10.000000, s: 0.333333) mcs: 20), ((beta: 10.000000, s: 0.444444) mcs: 20), ((beta: 10.000000, s: 0.555556) mcs: 20), ((beta: 10.000000, s: 0.666667) mcs: 20), ((beta: 10.000000, s: 0.777778) mcs: 20), ((beta: 10.000000, s: 0.888889) mcs: 20), ((beta: 10.000000, s: 1.000000) mcs: 20)]\n"
     ]
    }
   ],
   "source": [
    "import cxxjij.utility as U\n",
    "schedule_list = U.make_transverse_field_schedule_list(10, 20, 10)\n",
    "print(schedule_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "上の例では$\\beta=10$で固定しながら$s=0$から$s=1$まで、各パラメータで20モンテカルロステップ計算しながら10段階で$s$を変えていく設定例です。計200モンテカルロステップの計算を行います。\n",
    "$s$の変化については[Morita, Nishimori (2008)](https://aip.scitation.org/doi/10.1063/1.2995837)の手法を適用しています。\n",
    "\n",
    "#### Algorithmの実行\n",
    "\n",
    "続いて、Algorithmを実行します。前章と全く同じように書けます。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "import cxxjij.algorithm as A\n",
    "A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "前章と同じようにcallbackを使ってみましょう。横磁場イジング模型の場合は、システムとパラメータ (逆温度$\\beta$、$s$)を引数を持つ関数を作成すれば良いです。  \n",
    "例として、以下ではシステムのエネルギーの値を記録するcallbackを作っています。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "energies = []\n",
    "\n",
    "def callback_log_energy(system, t):\n",
    "    #graphは以前にGraphモジュールにて定義したオブジェクトです\n",
    "    #各trotterスライスの平均値から、古典スピンの0、1を決めます。\n",
    "    classical_spin = [-1 if np.mean(s)<0 else 1 for s in system.trotter_spins[:-1]] #最後のスピンは補助スピンのため、除く\n",
    "    energies.append(graph.calc_energy(classical_spin))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "このcallbackを用いて同じAlgorithmを実行します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "#スケジュールをもっと長く取ります (計20000モンテカルロステップ)\n",
    "schedule_list = U.make_transverse_field_schedule_list(10, 200, 100)\n",
    "A.Algorithm_SingleSpinFlip_run(mysystem, schedule_list, callback_log_energy)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "記録したシステムのエネルギーを、横軸をモンテカルロステップ、縦軸をエネルギーでプロットすると次のようになります。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": [
    "!pip install matplotlib"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAYEAAAEGCAYAAACD7ClEAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjQuMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8rg+JYAAAACXBIWXMAAAsTAAALEwEAmpwYAAA1U0lEQVR4nO3dd5wTdfoH8M+zuyxLW5Yufem9ulQREREQFDz1PPT0bCdnOctPT0WwgJ7C6XmW0zvFE8tZQQQ9UaRIR8AF6b0sAlIWpCxlF3b3+/tjJmGSTJJJMskkm8/79doXyWQy82w2zDPfLkopEBFRckpxOgAiInIOkwARURJjEiAiSmJMAkRESYxJgIgoiaU5HUAoatasqbKzs50Og4gooaxcufKwUqqW2WsJlQSys7ORm5vrdBhERAlFRHb7e43VQURESYxJgIgoiTEJEBElMSYBIqIkxiRARJTEmASIiJIYkwARURJjEnCIUgpTcvfg9Nlip0MhoiTGJOCQ95fm4ZHP1+KSF+c7HQoRJTEmAYcs3HYYAJBfUBTS+/b8ehpHT52NRkhkkVIKS7cfRmkpF2SixMck4BAJ830XvzAPXZ6dbWssFJrZGw/ixv8sx3tL85wOhShiTAJEIfrl2BkAwO4jpxyOhChyTAIOkXCLAkRENmIScAyzABE5j0nAIfFWEigpVfhuwwEoxcZOq/hJUVngWBIQkYYiMk9ENorIBhF5wKlYnOBEDvhqzS9YrPdK8vbe0jz86b8r8eXqX2IcVeKReMvgRBFwsiRQDOBhpVRbAD0B3CsibR2Mp8y7/5OfcNM7y01f2683dlrpsrrlQEFE3SNve3cFxn+7Kez3J6L/LtuN7FEzcKqIgwMpvjiWBJRS+5VSq/THBQA2AajvVDxkzfp9xzHolYX41/ztYR9j3pZ8vLVgZ8B9Nu0/gU37T4R9jnjz9kLt9z18MrRxIUTRFhdtAiKSDaALAJ/bVBEZKSK5IpKbn58f89iiJV5rFFSQmm5X98jVe45HNY4rXl2EK15dFNVzEFEcJAERqQxgKoAHlVI+t35KqYlKqRylVE6tWqbrJMe1o6fOInvUDEz7aa/Hdomz3kHRSEoFhecw/PXF2H6oIKLjnC0uxeo9x+wJykZsQ6eywNEkICLloCWAj5RSXzgZS7Ts0gcUvb/Uc53neC0JuGw7WBBwcruzJaU4UXgu4DEWbM3Hmr3H8fLsbRHF8tcZG3H1G0uwM/+kpf1LShWKiksiOmcg8f63s9O7S3bh3o9XOR0GRZGTvYMEwDsANiml/uFUHLHifdMYqwtJUXEJnpy+HsdOW5tvSCngXEkpLn95Ie7+0P9//oVb89Fx7Cy7wgxo3T6t6uno6cBJx+XOD3LR6omZ0QwpbLdMWoHhbyxxOgzLxv1vI2as3e90GBRFaQ6e+yIANwNYJyKr9W2jlVLfOBeS/fYd1erQ1zhUnTH9p33477LdKA7Sm8fY7bFE33fZziNRjS1avt98KKL3f7LiZ/RsWgNNalayKaLz8o6cBo6ctv24ROFysnfQYqWUKKU6KqU66z+OJYCVu4/ii1V7g+8Yoh0WqzC8FRWX4GxxacjvU0ph8o97cFLviui69oczCCzSKu9gpwy9m6m2f+G5Ehw8URheUBY8/sU6DPvnYgvRsFGAEp/jDcPx4tp/L8VDk9cA0KpDwrkAm/HXABysYbjVEzNx0d++9/v6icJz7gu90aJth/Ho1LV4avp6KKUwc/0BANYbMYPtFtZAKT9veXTqWmSPmoGl2w8je9QMzN54EA9PXoMPl+0O+PY7P8hFj+fnemwrtvFvBgAFAfrze8dzxauL8PgX6wIeL5naESixOFkdFJdufmc5FumjavMmDI34eH7/81u4KAQauNVx7CyIALvGD0Xe4VPo9/f5uD6nASbnaqWZ/JNFmLflEBZs9e1We/zMOVStUM5K+EHtyD+JSulpuKBqhs9rroRyqqgYB44XIiUFOHTi/O/0+Uot1r/P2gJAu7gDwNQgJbJFJqOeh72+BBv3n7DlbxYq15iG8dd0iPm5iSKVdCWB2RsP4rHP1/p93ewCE4lIbwCbj/4Gn6/ci8XbtLtlI6WAY6fPYv4WrQ7clQBcjpw0bwzuY1LCcCWcUOO97KUF6Dl+rulrriqo+Vvy0XP8XFz20gJcaVLNYkelysYoDCzzN7o3VpVA6/cdxwc/5MXobJSski4J3PlBLj7L3RPRMX45dsbvBeL4mXN4ZMqaoNMDWL3YFpcqTPh2E/63xnxOn87PmC8wIyIeFytj/XVBoW9s037ap+/n8SZbmZ0XCKGqKoJ4rPaOMmr39Hf4YUfkjeMlpcrd2B6KK/+5GE99uSHi8xMFknRJIFx5h0+5e8v0nvA9fvvmDwDgUy//7/k7MGXlXvR4fi7W7T3utzootLp1QUqAv5SVy8vavdZH+IZTfz1n48HQ36SLSjuEwfwth9D5mdl+J88LJDfvV994QjxGs9HfYK/eS4wo3jAJWNTv7/MxYuIy9/ON+0/gm3X70f7p7zDP0CXRVQVysqgYV73uv4eJ8UJyqEDr6VJSqvDFqr0+vWZSxPsdodt8wP+oXX+ljFD8Ua/PjyeutZhX7j4KAFj181G/+x45WWQ6+I39f6isYxKIwOhpWo+Q29770e8+/qowjDe3Hy77GQDwwQ95eGjyGkxasstj36LiUj0RmBv3v40+2xZuzce4r6xVJdz3yU+m8SooPPjpT1i0zdqcTQWGi2jhuRJsO2ixe6zFep6nvtzgTphWBFuL+aHPVuOVOVsBABf+dQ56Pe/btmEM7eCJQmSPmoGpq7SqM1fDNlEiY++gCJSU+F68IrlzdM0w+dcZntMsHz9zDj+aVEsEc+pseFMnGC9801f/gumrf0HehKEoPFeC/cf9V2uU6j00vRuw7bJx/wk8Nd2+OvIv9HaQBwe0BBD889p+SEtqrnmMCs/Z1yWVyClMAhGwcsE33vHP33IIx8+cw/DOvjNmt3tqZsCL0Fard9UhOHO2BOXTPAuDxgbkc15J7uIX5gXstjp/6yHUqlw+5DgCfY6//88yrN17zP28NEipobjE3gtzsAFh/gbhLd1xGDPXH8Azw9uHfM5DJwrdCYco2lgdBP9dAQGtS6lZX3vA/AKw/7hndYWxUfPWd3/EA5+u1rZ7Hijsu/ZwFRWXoM1TM31KHVsOFODJ6et99i8pVUEXnHng09W48T/mi9YEEui6vmT7EZ9kZOTdflIS4GB7fj2NXuPnuqfDDjU2s0O/s3iX70YAN769HB/8sNv0tWCGvb7E43Mc/MpCLN0eWdflBVvzkT1qBvIOn4roOFbYOWiPoi/pk8ChE4Xo+8I8v6/f+UEubpm0wv3cWNXh0QVTv0JYaWT9+chpx5coPKMnnSkrPbvLfrn6F0wxqesOpzoqGg4Ypou496NVaDr6G6zfd77n02c/+u/+O2XlXuw/XujuDmuF6288/ptN2HXY9+78Kxsa1b0d8JoSY/OBAjz55XqcKykNq6spoM0hBZg3ju89etq2aThW7j6Klk98a7kdiZyX9Elg4CsLceRU6H3IAc87w037C9wjXo22mPTK6fvivLhZTcBf331vxp5RTjJ2dZ2xTpvd0tiN00q/+pJSZX0uJaWw5UAB3lq4E09a7LMfzjxNLkcDfBdbjPkWV7y6MOxjm3nxu83o87d5PtNwhGvFLu1vsTjCkgvFTlIngbcW7MCxANMTB5pP39tPe45itklfeb93ioYs8Nr34S/VGK4xJlU+Tol0IrZA75676aBPNc4/Zm/F+0vzLB/7bzM3+33drI1iR/75KpczIVbzFQWpSjG2Da3dewzZo2YE7PoazBvzdoT93mg6dKIQN0xcFjApkj2SOgmM/9b/f24AQet0jRevRVtDu/NxelWqeJojPpqfxR3v52LOJt/kPG21tWqcdfuOB5yaev0+z+kqThUVY8A/Frifj/92k/dbAABTV+4NqcRgtuf8LVqVy/ebIps626ikVOGafy3xaQfbdrAArZ74Fou3Hca4/20IYwbY0ExcuBM/7DzCbrgxkLRJYN6W4P9xJgRJEsYugpsOhDZ3TbBeLmQfswbtonPW7tBD6aVz+mwx2j39nce2fX5GCj88ZQ0Gv7II/V+ab/n4LsfPaKXXUKoUrSacX0+dxaqfj+Hhyas9tn+yYg+Kiktx0zvL8e6SvLCnSA8Vp+uOvqRJAtmjZuDRz9e4n9/2rv8BXuEItY6fOeA8Jz4L4whqY7XfZK+G5VCme5hjckc+N0ApYsvBAuzMP+WeABCwdtEb+9UGvL/0/IXY+J7Hv1iLj5f/7POe6XrJ5/0fduN4gCrQT1f4vtdM0H4NEf5NOfV27CRNEgB8Z9m0U16Iq0WxG9152yJciN41yM6fYA3/bZ86f/f+6FT/M8wGc79h5HUobg3xhuRUUTGe/mqD+8JuTKKfrNiD0dPWeYytMFqz5xgGvKxVV3lX6Tw3YyNemr3VYhR+1snwc/E+V1KKy16aj7kmVXOB8GYp+pIqCcSTgiJr6+Umg0DjAKyI18bNUGw9GDgR7jQ0Ns/y6oBg9uk9oS8qVFxS6jOAzlU99oPX8qFvLzKOebD3VvzIybPYkX/KPdWKtwPHC3HzO8vPV3WxKBAzTAIOKWVBgAwGvxJ+18+DJwo9JjF0eWL6ejQf8y3+YXJ3X1xSir1HA5Veg61JHWqUgb0xbzsWbTuML1d7juFgQSD6OG2EQ7zvwii5RdLZ5otV+/DFqn3Y/twV7m1r9x53j6lwzXVkNP7bzX5HO7sEGyEeiL9fJ1j1jut1lgNihyUBojhx+GRRSFNaeGs+5lvL+wZLAICg23Nzgh7np5+PInvUDGzSV3YL9+Ltr2TBNoHoY0mAKE7k/DX4RTccS21YHc2b65o9c8MBANqYhTZ1MyM+rmJRIOZYEiAi2/i7cQ/WhuB62fv9Zl1mC8+V+KzoR+FjEiAiH8G63bp67/ibldTfNd9/kvB8hwQoCvT/+3y09xqUZ3TmbAm7YIeASYCIwvbdBq27qvf8SuFW5ftrA/jl2Bn35HS/HA8842mbp2baPtFeWZYUScDV95iIwvP5Sv9TdLt88EOe+/69uESZVtmEOtBYKe3/b+8J3+P6t36wFCvgOYkfBZYUScB7KgAiCs0JrynHzS7mxmm8Jy3Z5VFl47rD97cegnebgfF5p3GzQoqVQpMUScDVfY2I7PHo56FNrzElV7sRc03hMWvDAWSPmuEzUvrZrzd6rKJ27DSnko62pEgCR/lFIrLVCj8rzS3fZb7de/6mJ7/U1rMY+PJCvDBzM95dkud+7ZHP17pLGp5TWVA0JEUSmLeFS90R2e3dJb4X6EBrL7hs87r7/9f8xJ/7KZFxsBgRhWXc/zZa2q/j2O882hSsrITG+eNiJylKAkTkjLcW7PBpVAaAgyf8j0NQSnG6iBhiEiCiqHlt7jafbRMX7gz6Pn9VREXFoa3ZTMExCRBRTEXSh/+jZdZWPiPrmASIKGpOnbX3zr2QJQHbMQkQUVyxqzlAsWHBEkeTgIgMFpEtIrJdREY5GQsRxYf9QeYGsmrM9PW2HKescywJiEgqgDcAXAGgLYAbRKStU/EQUdny8XK2H1jhZEmgO4DtSqmdSqmzAD4FMNzBeIgozh0ydC192WTtZAqdk0mgPgDjzG579W0eRGSkiOSKSG5+Pkf+EiWz95bmuR+/atL9NJCi4hK8MW871xrwEvcNw0qpiUqpHKVUTq1atZwOh4gS1KTFeXjxuy1435BIyNkksA9AQ8PzBvo2IiLbnT6rjVw+c47dTI2cTAI/AmghIk1EJB3ACABfORgPEVHScSwJKKWKAfwZwHcANgGYrJTaEPhdRETnfbR8N25+ZznW7j0WdF8OGzDn6CyiSqlvAHzjZAxElLjGTNPGAuQXFGHmg30tvYcTlHriVNJEVOYNe30x1u497nQYcSnuewcREUWKCcA/lgSIKOGJCNbtPY6ConNOh5JwmASIqEy46vXFToeQkFgdREQJb9P+E5b35dKVnpgEiIiSGJMAEVESYxIgojJr7qaDyB41w2ObsD7IA5MAEZVZn6zY47OtoLDYgUjiF5MAEZVZJwp9u4y+uWCHx9KTxSWluPHtZVix69dYhhY3kjoJTL27t9MhEFEUWbmw7z9eiKU7juChyaujH1AcSuokUKl8qtMhEJEDFm8/jEKvKaX3Hj2DWRsOOBSRc5I6CaSwgYgoKd38zgo8+/VGn+0j/7vSgWicldRJoEG1Ck6HQEQO2Zl/ynS7SrI5p5M6CVRM56wZRORp0pI8p0OIKV4FiSgpicBnDAEA/LjrV9zRp4kDETkjqUsCdqtXNcPpEIiIQsIkYKPS5KpKJKIygEnARlUyWLtGlCjYOVDDJGCTCdd0QPPalZ0Og4gsEq42DIBJwDYjujdCkvUsI0po3oPFkhWTgI1KI8wCrS+oguGd69kUDREFkrv7qNMhxAUmARP//n3XsN4XacPwzAf74tURXYLuN25Yu8hORESWzd10EPO2HAIAnC0uxZPT1+PIySKHo7KPpSQgIveJSLVoB+OEV0d09ng+9qq2uKJD3bCOFWlJwKoODarG5DxEyWjmhgMY+9UG98jhO97PxW3v/ggA+Gbdfvx32W48980mJ0O0ldWSQB0AP4rIZBEZLGVoVYbhnet7PM+qmB72sUosFAXu7tfMdHvNyuXDPi8R2eu9pXmYnOu7FoHrRq8stf9ZSgJKqScAtADwDoBbAWwTkedFxPyKlqSslATKTPYkKuM2/OJ/8fqyNL+Q5TYBpf3WB/SfYgDVAHwuIi9EKTZHrR07EGueGhjSe4pLzL8Yb9wYXhsDETnH+zr/zbr9+PXUWWeCiSKrbQIPiMhKAC8AWAKgg1LqbgAXArg2ivHFnIL2l8/MKBfy4K8eTau7H3dtlOV+bK3yzJ47iys7hteeQUSelNf/yXs+WoW/zig7bQEuVksC1QFco5QapJSaopQ6BwBKqVIAV0Ytuhj5+M4elvc1u6Df3785AOC+/i3wwrUdAQDpaec/2rZ1M92P7bjUC4CrOpl3JX39xq7o1DDLhrMQJbdpq/aVqWoff6wmgVcBFIhIdcNPOQBQSiV8auzdrKbpdtefXwTY+fwQAMC1XRv47Hdp69oAgNQUwUUttGMNbncBqlYoBwDufwH7ehBdaChpePvoj9aTWjDP/aa9bceKpok3X+h0CFTGnDpbgtkbDzodRtRZTQKrAOQD2Apgm/44T0RWiUiZ/98nAFJSBGvHDsSEazr4vm4oHtTPqoCNzwzCLb2zDa+f37dR9Yphx5Fdw9p7K5e3bw6j3/dojIcub+n39bf/kGPbucJVPi0FLepUMX3NWCIjCtWJwmLT7WWpfGD1f8hsAEOUUjWVUjUAXAHgawD3APhXtIJzgvFGPUWAAW3qYNKt3QBo7QRpqcE/sorpafDXi/bG7o1Mt9evFvwCbzxmLHrpXthYGxpy/2UtTF/PmzAUF7cwL0UZVUpPxeZnB9sam9HrfhreK5dPw7R7eqMup/imKCkqLsGpIvNEkSisJoGeSqnvXE+UUrMA9FJKLQNQZju4iwj+c0sO+rWqbesxzUy6JfgdtWs5zIrpaRHVVTasHnhZzc3PDsZnI3ti6t29wzr+itGXeTwvl5aCjHKpfvd/Ymgb3NTTPDlaUS7V/DP9901d0a5eVSx69NKwj00UyPDXl6Dd098F3zGOWU0C+0XkMRFprP88CuCgiKQCKI1ifHFp+r0XeTz3d08+dlhbZGakWaqeqWFhsNj4azrgzZu6otUF5lUfRpXS/V90f9PFt10DAH6X0xAv/64TMsqlokfTGkHP4U/tTM8772D5KlCCsCpQuchf6e2RQa0iPi8lt80HCpwOIWJWk8CNABoAmA5gGoCG+rZUANdHJTKHWLnB7twwCxdkBq9i+E2XBlg7dpClKiQrKpdPw+D2WhfQTENjs5lVT12Ojc8MMn/Rzy/ZqEZFvwkiEsnQwwIAmtasFPVzfHiHfY3+FNzSHYdNt3+5+heP51ZmC4hXQa9O+t3+q0qp+5RSXZRSXfXH+Uqps0qp7TGIM2as3GUDwAd3dMegdnUwsm9TdKgf+7l8rvaa7sJb+bRUVExPwxND2+D927vbeu63/5CDJaP6A/A/BqJP8+BtBUauPPHs8NhNjmd3s4pxnEgw4U4TUr1S+NOaUOi+WLXP72vLdh5xP+77wrxYhBMVQZOAUqoEQGMRse3bJyIvishmEVkrItNEJMuuY5u51dBTxx/XXVyFANUoRi3rVMFbN+dg9JA2SEmJ7Gri6n7q7fO7emHzs4NRRa9O8mi0ThH0bVkr6LH/eHFTXOK1n5WqJ3/yJgzF5W3roH6W1q7gvTCHq/fTh4ZuqsHukTwuxmFcmUUk7At6lTB6UhkH5DWrZbz7DxzEmCFtTI8RCu8BTOScVT+fn4p637EzDkYSGav1FDsBLBGRJ0XkIddPBOedDaC9UqojtG6nj0dwrKBcVSf/+3MffH1fHwDAXwZ6dnvM1pOAHfXTofJOIuX1bo3t6lVFRrlUpPpp+PTeauVC+NoNXXBzz8Ye26be3RtVyqdhmJ8BaIGkp6Vg7FVt3c/T/MQKaA3GK8Zc5rO9fb2qEV3alFKmq0Rl1whePTPvkX6Y81DfkM734nWdMP3eizDj/j6oY6gWdDXc+9O1cZb7cVbFwNV54eISpxQqq0lgB7QuoSkAqhh+wqKUmqWUcvWrWgatvSHqOjSoivb1q2L1U5fj3kube7z2yojO+M8fctx3uE7yvphbrVLf8uwVfl977YYueOHajhjWqZ5P0rmwcTWsGzcIDcMcw3DrRU3cjwPlodqZGahdJQPXdNWqsmpXKY+1Ywd6jHC2q4ZmeOd6Hr/PZa19e3gJBDUrl0fz2qF9ldNSBZ0bZqFdvfPVgHf0aYI/9Goc4F3aGV28/6b+uuF6s/Jd6N0s/EZ9sq6sNHVZnUV0nFJqHIAXXY/153a4HcC3/l4UkZEikisiufn5+eGdweuvlVUx3aerZmZGOQxoWye844do+ejLsH6cn0ZbE65GJ38lApdAA6OGdaqH67s1tHzOcJk1WH82spfH839c3xlT7uqFr+/vg8yM0O+I77ok+OS1FdM974j7mSSBQK7PaYAFj/Qzfc34V2ipD1K7tmsDpJpUC/q7qfC+ftxt4XeyynvakOd/4zvAkcjF6gRyvURkI4DN+vNOIhJwkJiIzBGR9SY/ww37jIE2I+lH/o6jlJqolMpRSuXUqhW8Dtz0GLC/ETBUrS+o4p4+ok5mBiqXT0P9rAqmvYyeHd4emRlp7ot6+/ra3EPpXr2Mrs/RLurjr+mAKXd5XmhjzTU/0kUmU3C0rZfps61bdnXUrnL+d/d3V+XdngFog/iMwm0TCPSeKhnl0NhCddLoIW3w2cieaFsvExVMqhJrVPbTlKaUxxQkVtui/B5PJzj/vQCAEd0a4sYe4Y/BIP+m/eS/0TiRWK1AfAXAIABfAYBSao2IBKxIVUoNCPS6iNwKbfK5y1SU+xAq5fw8/jMf9P24XD1svP02pyF+a/iPPPEPOdh+6KRPe8XQjnUxtOPQiOLyF0OoujbOwsb9J1AnM7Kxg94X5vdv747Hv1iHT1b87Hcfc15fKZOvmB3fifS0FPeYChFBjUrpOGKYbtj4N6tdxfOzmXBtB0xdtdfyuRY9einqVg1eXen6vdJSBONNpjkhe2w/dNLpEGwRynoC3svslIR7UhEZDOBRAMOUUqfDPY5VCiom0yxES2ZGOXRtFJ3VPeOhDQQAKpfXLpYZab53xJkVtHuVprWs98P3vuab3WV0bez/M7X72/LidR092igUgHIBxo+sHzfIp8rQrM3ms5E9/R6jfrUKCf29p9iwWhLYIyK9ASh99tAHAEQye+jr0KabmK1/SZcppe6K4HgBxUNJoKy7qWdjTF25L+x2lYcub4VqldJxdZf6+MfsrR5d7v5vQEvUqlweJ4uK8cqcbT49gUJJZDf3bIxnr26PU0XFqBTmRHuBLqzeyWbsVe0w9qsNPlN/Byv7Whll/uJ1HdGjaQ30aFIdy3f9CkDrAFAvqwKa166MJ69sG+QIRNaTwF3QppOuD2AfgFkA7g33pEqp5sH3sk88tAmUda0vyMSmCCaJq5Ceinv6aV+Laff2xoc/7MZv9DrzjHKp+OPFTfHZj1qVUJ3M8hjRrSEm5+7B4sf6o15WBez5NXCBsoXeA6hjA61HT7AEcHWXwIPxrGpbLxOTo9ReY6wydHHNczXnoUs8tk/+Uy8s2paPf35fpsZ2kg0sJQGl1GEAv49yLFGjlQSYBRJF7SoZeGig77w+1+c0RNUK5TCw7QVISRFM0BfwMeN9p92rWQ0seKSfpam8//OHHLSP8ijwQA28oa7h4LrB+SDAyPDuTaqje5PqYSWB12/sgj9//FPI76PEYCkJiEgtAHcCyDa+Ryl1e3TCspcC64O8fTayp+UeKfFCRNxzJ3mrl1UBQzpoC/l8ssK7+UpjpbePpTjCfM3oD72y/b5m1svIikAD9aziBT98328+iP6tY9PN3E5WG4a/BFAVwBwAMww/iYE5wEePpjXQsUGW02HYJjVF8K/fX4hO+u9k1/QKrhHmVlk9q9mYAjNm6zDMffgS2+eDcmkShUnwFj16qceo8rLq9vdynQ4hLFbbBCoqpR6LaiRRxDaB5GH33zka1UKBpvkGPKuyzKYxaVarMprVqux+/sTQthj1xVp0DnNt6cvb1onqMooNq1dE3TjphUa+rJYEvhYR81nOEoC/uWWIzFidSTYcU+/ujbkP9wvpPSP7NsWlrfwPlGxfvyq+vu9in1HSVv3fgPPzaCmltYlQ8rD6rXkAwOMichbAOWi1K0op5TsUNE6xJJBcwh1+mDch+OA7q9+lwe0u8Nl2odfYhDdu7IqTRec8tnmHPtow+6hdMsqloPCcth6U94juUKbEBrRBacUJPJ++nQ4VFHqMhE8EVksCVQHcCmC8fuFvB+DyaAVlt7Iy0RMF5yrxOfUnf8awHoJrkFsgQzvWxe+6xX5ah4tbeJYsBrTRupaaVT8FK0VbWeMgWe7Buj831+kQQmY1CbwBoCeAG/TnBdAGfCUEheT5Eia9KP6hb+ge/GJ9Zcd6+Nu1kU3V4JorKpq8P6aXru+MV0d0RvPalU33DyTUCerMFmGy2lCeCL7ffBD//SEPQ15dlBDrDFitDuqhlOoqIj8BgFLqqJ2LzESbUoFHeZI11UKcA/+jP/aIeTVc10ZZAIAr2vtWxQRycYuaWLTNfCnBbtnV0LZuJp6+qh2eurJt1L5Lq5+6HMfPnLOtK2soqlYoh+H6anX+fr/yaSkoKvZdUryyhTUMjMdsXKMi1u077vF6ZkYajp4+5/22hGTsJfT+0ryoVOfZyWoSOKcvM6kA97iBhFlgXkGxJBChtWMHIi3Eu7WLQlxi0g7Na1exVK/vbeLNOcgvKDJ9bcpdvd2PrYyt6KNXtYywUHIwyqqYjqyK8XlvNfXu3qiXlYFe47/3ec3KUpn8/xe/rCaB16AtMF9bRJ4DcB2AJ6IWlc0U64MiFs68/4mkQnoqGtUIb1Edb/WzKoSViGLpmeHtMSuEbqHeDdqANs33Y4Nbh1WF5I0ldedYnTbiIxFZCeAyaJfTq5VSkUwgF3P8ihGdd0FV6z1Y/F2f37zpwoQbdR5riXDdsdyxWCm1GfqiMolGqcSeSpooHoWSAIz//cz+L0Z5SREKwPJ6AomMI4aJYiejXAqyvarWkvX/31sLdyJ7lPkMO2fOluB4HDSGJ0cS4NxBRLZY+UTABQMBAJNu7eZeAMfVS6tHE231tQbVKpiu/5CMJfXLX16ATs/McjqMJEkCCb6yGJGTjP9zagTpCVSjUjo61K+K8voKca51GSqVT0PehKFY/Fj/gF2Nn06CieZc9h6NjzEESZEE0lJSfBZpJ0p2n47siX/e0MVne7CewDkBluVc+eTlqGLoSWZ2qNsuagIAGDesnc9rw7xWYKPoC2/GqQQzdlg7jDX5whElmxWjL8OJwmIAQM+mNUz3qZiehheu7YgPl+/G2r3HfV7/6M4eKDwb/jCh9LQUdxfap7/aAADonl0dMzccQPkw11Kg8CVFEiAiTe3MDNS2MCvF9d0aYu2+Y1i79zgqei3FWT4t1V3d49KwegUcPGE+2M6KV0Z0xu4jpy2trUz24idORKbGDGmLDvWrom+L4CO/5//l0ojOlVEu1T2F96ODW+H9pXmolJ6GnYdPRXTcePHLsTOoF6drKrCinIhMVUhPxe+6NbLUqSI1RTwmgXv6qrYY1K4O+rb0vw6CP/f0a47lowfg1ouyfV6765JmIR8vHoz/Nn6HWDEJEJHtGlaviLduzjGdmjoSl4SRVOLB/9b84nQIfjEJEJFjPri9O0Z0axhwH+OFv1cz88bsRJA9agaOnAy/3SRa2CZARI7p27JW0CqjRtXtmdgvHmw5WIDeFmZdjSWWBIiIYiQep0hiEiCiuOSaamKQyVrNierj5T+j6eMzUFwSP8uxsDqIiOJSqwvCWyAons1Ytx8AUGiyQptTWBIgIoqxeJrJjEmAiBLS6CGtnQ4hbPE0nyWTABElpJF9E3PgWLxhEiCihLV27EDT7Wazo8YTiaMKISYBIkpYmYZpq1vVqeJ+fFWneriyY10nQrIk70j8zInEJEBEca92lfK4r3/zgPt8OrInAKBqBS0xxGGXfLc1e445HYIbu4gSUdxbMSb4spY+ja1xnAXeWbzL6RDcWBIgooTWtVEWAKBKRjnkNK6GV0Z0djQeK7YdOul0CG6OlgRE5GEAfwdQSyl12MlYiCgxTf5TLxSXKqSmCD6/u7d7u4rnokAccSwJiEhDAAMB/OxUDESU+NJSU5BmMmN1PM7TE4+crA56GcCjiOuaOyJKVN5JIG/CUAxoU9uZYAJQDmcrR5KAiAwHsE8ptcbCviNFJFdEcvPz82MQHRFR7MzbcsjR80etOkhE5gAwm/5vDIDR0KqCglJKTQQwEQBycnJYaiBKcnf3a4ZZGw4E3c+8TSB+Bmm55Bc4u9BM1JKAUsq0T5eIdADQBMAafe3SBgBWiUh3pVTwvywRJbXHBrfGY4ODzxtkVstSt2pGFCJKbDGvDlJKrVNK1VZKZSulsgHsBdCVCYCI7GRWDhgztE3M4wjG6SkkOE6AiJKG98L3T17Z1qFIznO6K6vjSUAvEXCMABHZykqnm/6t46+3UKw5ngSIiKIjeBZwzTPkJKfHMzAJEFGZ5O/i+uCAFu7HGeWcvwQ6vcCM858AEVEMldeHF9/cszEqpnMOTSYBIiqTbr0o23T7bRdl477+zX16CvVqWiMGUflidRARURRc3KIWHhnUCgAwuN35casZ5VLx8MBWPj2FGteoGNP44gWTABGVeU1qVQq6zwOGtoJYqlTe2SopJgEiKvOsVLnUrVoh+oGYmLJyryPndWESIKIyK9SeN92bVI9OIAEs3OrsxJhMAkRU5lkdlTv5T70wqF2doPv1bVkr0pDiBpMAEZVZ4czL89bNOcibMDTgPv1bMQkQESUtcXqEl42YBIio7LO5L34ZygFMAkRUdpWli3W0MAkQEYWoLOUWJgEiKvNsn5mhDBUxmASIqMyK1qU6Wsc9WVSMgsJzUTq6OSYBIqIAmtY0n3LC3/ZgUlPMU8jmAyfQ/unv0GHsrLCOGy4mASIq81QEU3V+/5d+ptvnPHRJWMdLNalKWr7zCO7/5KewjhcpJgEiKrOiVXUfyXHPlpT6bPvn99ux9eDJCCIKH5MAEVEQX957Ed64sav7eSgjkYd0uMBnW+eGWR7PSx1cVIBJgIjKrP6ttXmArupUL6LjdGqYhaEd67qfd22cZbk0UL1Sus+2j+/sgXmGaiYnkwDXViOiMqt57cpB5wEKR+sLMi23M5iVGiqmp6Fe1vntG385YVtsoWJJgIjIxK29s9GxQdWIj+OvxGBsID5RWOzx2nMzNuL02WLvt0QFkwARkYmxw9rhqz/38dles7JWvSMilkoZ/hazTwlQn/T2ol349/wdFiONDKuDiIgs+vq+Prigaobl/T+/qxfW7ztu+lqKn/ECLkXFvr2IooFJgIiS2ru3dcPaPeYXam/t64dWPZSTXR3r/CSBYEpKY9NYzCRAREnt0la1cWmr2jE7X5dGWZb2i1WPIbYJEBHF0M09G1vaL1a9RpkEiIgisPnZwaiUnoqFj1yKYYbxCFPu6mW6f6AGYaNYVQcxCRARRSCjXCo2PDMYjWpUxGs3dAEAtKmbiW7Z1U33tzrILFbVQWwTICKy0df39UHDahUjPs6ZcyU2RBMckwARkY2C9SDKrFDO0nHyC4rsCCcoJgEiohgY1qkemtaqhEta1LK0v791B+zGNgEiohioXikdDw5oGXSQmMv8Lfn4+cjpKEfFkgARUVRF0r7b98V5qKePUH7p+s7o1ayGTVGdxyRARBQn+rWqhflb8t3POzXMQsvalQEA1SpZa0sIlWNJQETuA3AvgBIAM5RSjzoVCxFRPCiX6llDP/qK1ujR1P67fyNHkoCIXApgOIBOSqkiEYndmG0iojjVqLpn19LTMegm6lTD8N0AJiiligBAKXXIoTiIiOJGw2oVsPP5Ie7nvaJcCgCcSwItAVwsIstFZIGIdPO3o4iMFJFcEcnNz8/3txsRUVz6TZf66NIoC3f2bRp0XwXPKaYzyqVGMTJN1KqDRGQOAN8VloEx+nmrA+gJoBuAySLSVJms16aUmghgIgDk5OQ4txAnEVEYqlVKx7R7LrK0b4UYXPS9RS0JKKUG+HtNRO4G8IV+0V8hIqUAagLgrT4RJaWmtSrhugsbAADSUgQV0mOTEJzqHTQdwKUA5olISwDpAA47FAsRkeOuz2mINL130Fd/7hOzEcNOJYFJACaJyHoAZwHcYlYVRESULDbtP+F+3LZeZszO60gSUEqdBXCTE+cmIopHsbrz98a5g4iI4kDFGLUBeGMSICJyyJ0XN3E/dqJnEMAkQETkmDFD26JvS21q6ayK6Y7EwCRAROSgt266ECP7NsUdfZoE3zkKOIsoEZGDKqSnYvSQNo6dnyUBIqIkxiRARJTEmASIiJIYkwARURJjEiAiSmJMAkRESYxJgIgoiTEJEBElMUmkGZxFJB/A7jDfXhPxuWYB4woN4woN4wpNvMYFRBZbY6VULbMXEioJREJEcpVSOU7H4Y1xhYZxhYZxhSZe4wKiFxurg4iIkhiTABFREkumJDDR6QD8YFyhYVyhYVyhide4gCjFljRtAkRE5CuZSgJEROSFSYCIKIklRRIQkcEiskVEtovIqCifq6GIzBORjSKyQUQe0LePFZF9IrJa/xlieM/jemxbRGRQNOMWkTwRWafHkKtvqy4is0Vkm/5vNX27iMhr+vnXikhXw3Fu0fffJiK3RBBPK8NnslpETojIg059XiIySUQOich6wzbbPh8RuVD//Lfr75UI4npRRDbr554mIln69mwROWP47N4Mdn5/v2OYcdn2txORJiKyXN/+mYhYWoPRT1yfGWLKE5HVDnxe/q4Pzn3HlFJl+gdAKoAdAJoCSAewBkDbKJ6vLoCu+uMqALYCaAtgLIC/mOzfVo+pPIAmeqyp0YobQB6Aml7bXgAwSn88CsDf9MdDAHwLQAD0BLBc314dwE7932r642o2/a0OAGjs1OcFoC+ArgDWR+PzAbBC31f0914RQVwDAaTpj/9miCvbuJ/XcUzP7+93DDMu2/52ACYDGKE/fhPA3eHG5fX6SwCecuDz8nd9cOw7lgwlge4AtiuldiqlzgL4FMDwaJ1MKbVfKbVKf1wAYBOA+gHeMhzAp0qpIqXULgDb9ZhjGfdwAO/rj98HcLVh+wdKswxAlojUBTAIwGyl1K9KqaMAZgMYbEMclwHYoZQKNCo8qp+XUmohgF9Nzhnx56O/lqmUWqa0/60fGI4VclxKqVlKqWL96TIADQIdI8j5/f2OIccVQEh/O/0Otj+Az+2MSz/u9QA+CXSMKH1e/q4Pjn3HkiEJ1Aewx/B8LwJflG0jItkAugBYrm/6s16km2QoPvqLL1pxKwCzRGSliIzUt9VRSu3XHx8AUMeh2EbA8z9mPHxegH2fT339cTRivB3aXZ9LExH5SUQWiMjFhnj9nd/f7xguO/52NQAcMyQ6uz6viwEcVEptM2yL+efldX1w7DuWDEnAESJSGcBUAA8qpU4A+DeAZgA6A9gPrTjqhD5Kqa4ArgBwr4j0Nb6o3z3EvN+wXtc7DMAUfVO8fF4enPp8AhGRMQCKAXykb9oPoJFSqguAhwB8LCKZVo9nw+8Yl387gxvgebMR88/L5PoQ0fEikQxJYB+AhobnDfRtUSMi5aD9gT9SSn0BAEqpg0qpEqVUKYC3oRWBA8UXlbiVUvv0fw8BmKbHcVAvRrqKwIcciO0KAKuUUgf1+OLi89LZ9fnsg2eVTcQxisitAK4E8Hv94gG9uuWI/ngltPr2lkHO7+93DJmNf7sj0Ko/0kziDYt+rGsAfGaIN6afl9n1IcDxov8ds9KYkcg/ANKgNZo0wflGp3ZRPJ9Aq4d7xWt7XcPj/4NWNwoA7eDZWLYTWkOZ7XEDqASgiuHxUmh1+S/Cs1HqBf3xUHg2Sq1Q5xuldkFrkKqmP64eYWyfArgtHj4veDUU2vn5wLfRbkgEcQ0GsBFALa/9agFI1R83hXYRCHh+f79jmHHZ9reDVjI0NgzfE25chs9sgVOfF/xfHxz7jkXlQhhvP9Ba2LdCy/BjonyuPtCKcmsBrNZ/hgD4L4B1+vavvP6jjNFj2wJDS77dcetf8DX6zwbXMaHVvc4FsA3AHMOXSQC8oZ9/HYAcw7Fuh9awtx2Gi3eYcVWCdtdX1bDNkc8LWjXBfgDnoNWn3mHn5wMgB8B6/T2vQx+1H2Zc26HVC7u+Z2/q+16r/31XA1gF4Kpg5/f3O4YZl21/O/07u0L/XacAKB9uXPr29wDc5bVvLD8vf9cHx75jnDaCiCiJJUObABER+cEkQESUxJgEiIiSGJMAEVESYxIgIkpiTAIU90REiciHhudpIpIvIl+HebwsEbknjPdVFpG3RGSHPu3GfBHpEeIx5otIRIuFi8jVItI2kmMQuTAJUCI4BaC9iFTQn1+OyEaOZgEIOQkA+A+0SclaKKUuBHAbgJpW3ywiqWGc08zV0GaeJIoYkwAlim+gjZ4EvOZ+0edin65PWLZMRDrq28fqE5jNF5GdInK//pYJAJrpc8e/qO/7iIj8qB9jnPfJRaQZgB4AnlDadAhQSu1SSs3QX5+ulw42GCbmg4icFJGXRGQNgF5ex7xBn/d9vYj8zeyXFpEJos09v1ZE/i4ivaHNsfSiHn8z/Wemfv5FItJaf+97IvKmiOSKyFYRuTLUD52SQCQjPfnDn1j8ADgJoCO0KYUzoI2y7Afga/31fwJ4Wn/cH8Bq/fFYaFNjlId2x34EQDn4TnMwENoi3gLtxuhrAH29YhgGYFqAGF0jPCtAG61ZQ3+uAFxv2G8+tBGd9QD8DG3KgjQA3wO42uuYNaCNrHUN6szS/30PwHWG/eZCK50AWqL63rDfTP13agFt5GyG039P/sTXj2tiJqK4ppRaq0+9ewO0UoFRH2hD/6GU+l5EahhmgZyhlCoCUCQih2A+5e9A/ecn/XllaBfNhSGEeL+I/EZ/3FB//xEAJdAmC/PWDcB8pVQ+AIjIR9AWQplu2Oc4gEIA7+jtHz5tIPpslL0BTJHzC0iVN+wyWWkll20ishNAa2hJlAgAmAQooXwF4O/QSgE1LL6nyPC4BObfeQEwXin1VoDjbADQSURSlVIlHm8W6QdgAIBeSqnTIjIfWokFAAq997dKKVUsIt2hLbZzHYA/QyvpGKVAm3O/s7/DBHlOSY5tApRIJgEYp5Ra57V9EYDfA+4L8mHlNUe7lwJoS/u5fAfgdv2uGiJSX0RqG9+glNoBIBfAOBH3OrPZIjIUQFUAR/UE0BraDI7BrABwiYjU1BuMbwCwwLiDHk9VpdQ30Gbj7OQdv/577hKR3+rvERHpZDjMb0UkRW/TaAqteonIjSUBShhKqb0AXjN5aSyASSKyFsBpALcEOc4REVki2iLk3yqlHhGRNgB+0K/vJwHcBN854v8IbYGU7SJyBsBhAI9AmxHyLhHZBO0iu8zC77JftAXV50EricxQSn3ptVsVAF+KSIa+z0P69k8BvK03dF8HLQH+W0SegNbm8Sm0mWIBrd1hBYBMaLNnFgaLjZILZxElKqNE5D1ojeefB9uXkherg4iIkhhLAkRESYwlASKiJMYkQESUxJgEiIiSGJMAEVESYxIgIkpi/w+A+IZCm7vI/AAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "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": [
    "## 結果の取得 -Result-\n",
    "\n",
    "`result.get_solutions`で計算結果である古典スピンを取得します。この関数は最適化問題を解く観点にフォーカスを当てているため、trotterスライスの中でもっともエネルギーが低いスピン列を返します。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "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": [
    "## C++ core interface"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "```cpp\n",
    "#include <graph/all.hpp>\n",
    "#include <system/all.hpp>\n",
    "#include <updater/all.hpp>\n",
    "#include <algorithm/all.hpp>\n",
    "#include <result/all.hpp>\n",
    "#include <utility/schedule_list.hpp>\n",
    "#include <utility/random.hpp>\n",
    "#include <random>\n",
    "\n",
    "#include <iostream>\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<double>(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; i++){\n",
    "        for(std::size_t j=0; j<N; j++){\n",
    "            dense.J(i, j) = (i == j) ? 0 : gauss(rand_engine)/N;\n",
    "        }\n",
    "    }\n",
    "\n",
    "    //set local fields\n",
    "    for(std::size_t i=0; i<N; i++){\n",
    "        dense.h(i) = gauss(rand_engine);\n",
    "    }\n",
    "\n",
    "    //create transverse Ising system\n",
    "    auto system = system::make_transverse_ising(dense.gen_spin(rand_engine), dense, 1.0, 4);\n",
    "\n",
    "    //generate schedule list\n",
    "    auto schedule_list = utility::make_transverse_field_schedule_list(10, 20, 10);\n",
    "\n",
    "    //do annealing (updater: SingleSpinFlip)\n",
    "    algorithm::Algorithm<updater::SingleSpinFlip>::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",
    "```"
   ]
  }
 ],
 "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": 2
}