notebooks/5.1_variational_quantum_eigensolver.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "## 5-1. Variational Quantum Eigensolver（VQE）アルゴリズム\n",
    "\n",
    "まず、物質科学と量子化学への応用が期待されているアルゴリズム、VQE (Variational Quantum Eigensolver : 変分量子固有値ソルバー) アルゴリズムを紹介する。このアルゴリズムは物質の基底エネルギーの値を求めるのに用いられる。\n",
    "\n",
    "### 背景\n",
    "\n",
    "分子や物質の性質は、ほとんどの場合その中に存在する電子の動きによって決まっていると考えられている。したがって、電子を支配する方程式であるシュレディンガー方程式([4-1節](4.1_quantum_simulation.ipynb)も参照)\n",
    "\n",
    "$$ H|\\psi\\rangle = E|\\psi\\rangle$$\n",
    "\n",
    "を解けば、分子や物質の性質を計算によって明らかにすることができる。ここで $H$ はハミルトニアンと呼ばれる演算子 (行列) であり、分子の形など、系の詳細によって決まる。 $H|\\psi\\rangle = E|\\psi\\rangle$ の形からもわかるように、シュレディンガー方程式を解くということは、ハミルトニアン $H$ の固有値問題を解き、固有値 $E_i$ と対応する固有ベクトル (固有状態とも呼ばれる) $|\\phi_i\\rangle$ を求めることと同値である。このとき固有値 $E_i$ は、固有状態 $|\\phi_i\\rangle$ の持つエネルギーとなる。\n",
    "\n",
    "この問題を解くことを、化学分野では**量子化学計算**と呼び、現在も問題の解法について活発な研究が進められている。というのも、この問題は電子の数に対して指数的に難しくなるため、分子のサイズや結晶の単位格子が大きくなれば、厳密に解くことは実質不可能になる。そこで様々な近似解法が研究されているのだ。この問題の難しさは、$H$ が量子状態に作用する行列であって、粒子数が多くなればその次元が指数的に大きくなっていくところを原因としている。\n",
    "\n",
    "### 変分法\n",
    "極限的な環境を考えない限り、電子の状態は通常一番エネルギーの低い状態、すなわち基底状態に落ちていることがほとんどである。そこで固有状態の中でも、特に基底状態が興味を持たれる。\n",
    "\n",
    "非常に大きな次元を持つハミルトニアンの基底状態を求めるのに有力な手法として、**変分法**がある。変分法では、任意の状態 $|\\psi\\rangle$ についてそのエネルギー期待値が必ず基底エネルギー $E_0$ よりも高くなる、つまり\n",
    "\n",
    "$$\\langle \\psi|H|\\psi\\rangle \\geq E_0$$\n",
    "\n",
    "となることを使う。(このことは[変分原理](https://ja.wikipedia.org/wiki/変分原理#リッツの変分原理)とも呼ばれる。) このことから、ランダムに状態 $\\{|\\psi_i\\rangle\\}$ をたくさん持ってきて、その中で一番エネルギーが低い状態を見つければ、それは $\\{|\\psi_i\\rangle\\}$ の中で最も基底状態に近い状態であると考えられるだろう。\n",
    "\n",
    "実際には、ランダムに状態を持ってきていたのでは、基底状態に近い状態が得られる確率は系のサイズに対して指数的に小さくなってしまう。そこで普通は物理的・化学的直観や経験をもとにパラメータ付きの量子状態 $|\\psi(\\theta)\\rangle$ ($\\theta$ はパラメータ) を構成し、\n",
    "\n",
    "$$\\langle \\psi(\\theta)|H|\\psi(\\theta)\\rangle$$\n",
    "\n",
    "を最小化するような $\\theta$ を見つけるというアプローチがとられる。古典コンピュータ上で計算をする都合上、これまでは量子状態 $|\\psi(\\theta)\\rangle$ は、古典計算機でも効率的に記述できるものの中から選ぶ必要があった。\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### VQEとは\n",
    "\n",
    "VQE とは、**変分法に対して、量子コンピュータで効率的に記述できる量子状態を用いて基底状態を探索するアルゴリズムである。**\n",
    "\n",
    "アルゴリズムは以下のようになる。\n",
    "\n",
    "1. 量子コンピュータ上で量子状態$|\\psi(\\theta)\\rangle$を生成する。\n",
    "\n",
    "2. $\\langle H(\\theta)\\rangle = \\langle \\psi(\\theta)|H|\\psi(\\theta)\\rangle$ を測定する。  \n",
    "\n",
    "3. 測定結果をもとに、古典コンピュータによって $\\langle\\psi(\\theta)|H|\\psi(\\theta)\\rangle$ が小さくなるような $\\theta$ を決定する。\n",
    "\n",
    "これを $\\langle\\psi(\\theta)|H|\\psi(\\theta)\\rangle$ が収束するまで繰り返すことで、近似的な基底状態を求める。  \n",
    "\n",
    "![VQE-concept](figs/5/VQE_concept.png)\n",
    "\n",
    "(図の引用：参考文献[1])\n",
    "\n",
    "\n",
    "### 実装例\n",
    "\n",
    "以下ではnumpyを使ったVQEの実装例を示す。（もちろん、Qulacsでも簡単に実装できるので、余裕のある読者はチャレンジしてみてほしい）   \n",
    "ここではH-He$^+$(水素化ヘリウムイオン)の基底エネルギーを求める。使用する量子ビットは2個で、参考文献[1]に従って実装していく。\n",
    "\n",
    "### 量子ゲートの準備"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "import numpy as np\n",
    "nqubits = 2\n",
    "#パウリ演算子を準備する。\n",
    "pI = np.array([[1+0.0j,0+0.0j],[0+0.0j,1+0.0j]])\n",
    "pX = np.array([[0+0.0j,1+0.0j],[1+0.0j,0+0.0j]])\n",
    "pZ = np.array([[1+0.0j,0+0.0j],[0+0.0j,-1+0.0j]])\n",
    "pY = np.array([[0+0.0j,-1.0j],[0.0+1.0j,0.0+0.0j]])\n",
    "pHad = (pX+pZ)/np.sqrt(2)\n",
    "pP0 = (pI+pZ)/2\n",
    "pP1 = (pI-pZ)/2"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "#パウリ演算子を1量子ゲートに変換する。\n",
    "X=[1]*(nqubits)\n",
    "Y=[1]*(nqubits)\n",
    "Z=[1]*(nqubits)\n",
    "H=[1]*(nqubits)\n",
    "P0=[1]*(nqubits)\n",
    "P1=[1]*(nqubits)\n",
    "\n",
    "for i in range(nqubits):\n",
    "    for j in range(nqubits):\n",
    "        if(i != j):\n",
    "            X[i] = np.kron(pI,X[i])\n",
    "            Y[i] = np.kron(pI,Y[i])\n",
    "            Z[i] = np.kron(pI,Z[i])\n",
    "            H[i] = np.kron(pI,H[i])\n",
    "            P0[i] = np.kron(pI,P0[i])\n",
    "            P1[i] = np.kron(pI,P1[i])\n",
    "        else:\n",
    "            X[i] = np.kron(pX,X[i])\n",
    "            Y[i] = np.kron(pY,Y[i])\n",
    "            Z[i] = np.kron(pZ,Z[i])\n",
    "            H[i] = np.kron(pHad,H[i])\n",
    "            P0[i] = np.kron(pP0,P0[i])\n",
    "            P1[i] = np.kron(pP1,P1[i])\n",
    "Ide = np.eye(2**nqubits)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "#2量子ゲートを準備する。\n",
    "CZ = [[0 for i in range(nqubits)] for j in range(nqubits)]\n",
    "CX = [[0 for i in range(nqubits)] for j in range(nqubits)]\n",
    "\n",
    "for i in range(nqubits):\n",
    "    for j in range(nqubits):\n",
    "        CZ[i][j]= (P0[i]+np.dot(P1[i],Z[j]))\n",
    "        CX[i][j]= (P0[i]+np.dot(P1[i],X[j]))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "#変分量子ゲート(X,Y,Zに関する回転の角度を指定できるゲート)を準備する。\n",
    "from scipy.linalg import expm\n",
    "def RX(target,angle):\n",
    "    return expm(-0.5*angle*1.j*X[target])\n",
    "def RY(target,angle):\n",
    "    return expm(-0.5*angle*1.j*Y[target])\n",
    "def RZ(target,angle):\n",
    "    return expm(-0.5*angle*1.j*Z[target])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "#初期状態|0000・・・0>を準備する。\n",
    "def StateZeros(nqubits):\n",
    "    State = np.zeros(2**nqubits)\n",
    "    State[0]=1\n",
    "    return State"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### ハミルトニアンを準備する\n",
    "\n",
    "参考文献[1]の [Supplementary Information](https://static-content.springer.com/esm/art%3A10.1038%2Fncomms5213/MediaObjects/41467_2014_BFncomms5213_MOESM1050_ESM.pdf) の表から、H-He間の距離が$0.9$オングストロームの時のハミルトニアンの係数を読み取り、定義する。このハミルトニアンの最小エネルギー固有状態を求めれば、様々なH-He$^+$分子の性質を知ることができる。\n",
    "\n",
    "※ このハミルトニアンは、電子-原子核間のクーロン相互作用および電子同士のクーロン相互作用の大きさから導かれている。詳細は第6章の量子化学計算に関する項で学ぶことになる。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "M = (-3.8505 * Ide  - 0.2288 * X[1] - 1.0466 * Z[1] - 0.2288 * X[0] +  0.2613 * np.dot(X[0],X[1]) + \\\n",
    "     0.2288 *np.dot(X[0],Z[1]) - 1.0466*Z[0] + 0.2288* np.dot(Z[0],X[1]) +  0.2356 * np.dot(Z[0],Z[1]) )/2"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 量子回路を準備する\n",
    "\n",
    "論文と全く同じ形式の変分量子回路を以下のように実装する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "n_param = 6\n",
    "def TwoQubitPQC(phi):\n",
    "    state = StateZeros(2)\n",
    "    state = np.dot(RX(0,phi[0]),state)\n",
    "    state = np.dot(RZ(0,phi[1]),state)\n",
    "    state = np.dot(RX(1,phi[2]),state)\n",
    "    state = np.dot(RZ(1,phi[3]),state)\n",
    "    state = np.dot(CX[1][0],state)\n",
    "    state = np.dot(RZ(1,phi[4]),state)\n",
    "    state = np.dot(RX(1,phi[5]),state)\n",
    "    return state"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 量子状態のエネルギー期待値を測定する\n",
    "\n",
    "変分量子回路によって出力される状態のエネルギー期待値を算出する関数を以下のように定義する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "def ExpectVal(Operator,State):\n",
    "    BraState = np.conjugate(State.T) #列ベクトルを行ベクトルへ変換\n",
    "    tmp = np.dot(BraState,np.dot(Operator,State)) #行列を列ベクトルと行ベクトルではさむ\n",
    "    return np.real(tmp) #要素の実部を取り出す"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### エネルギー期待値の最小化\n",
    "エネルギー期待値の最小化を、`scipy.optimize.minimize` に実装されている Powell 法によって行う。Powell 法は勾配情報を使わない最適化手法の一つである。パラメータの初期値はランダムに指定する。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZQAAAEKCAYAAAA1qaOTAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl8VfWd//HXJxv7TkBWgZBgK0XEiIgbylLrqNjWtmptsXZ+aq2KSH9jO/39pjO/6XQ6WqCodcGqZTpUbV1a60pAcUcNFAEFwqpsQgDZ1ySf3x/3xIZ4kxzIvffcm7yfj8d55J5zv+ee9/FiPvme5XvM3REREWmsrKgDiIhI06CCIiIiCaGCIiIiCaGCIiIiCaGCIiIiCaGCIiIiCaGCIiIiCaGCIiIiCaGCIiIiCZETdYBU6tq1q/fr1y/qGCIiGWXBggXb3D2/oXbNqqD069eP0tLSqGOIiGQUM/soTDsd8hIRkYRQQRERkYRQQRERkYRQQRERkYRQQRERkYRQQRERkYRQQRERkYRQQQnh7dXbuXfeqqhjiIikNRWUEF5evoVfvbSCVVv3RB1FRCRtqaCEcMN5BbTMzWbanJVRRxERSVsqKCF0aduCa8/qz3OLN/Phpt1RxxERSUsqKCH9r3MG0K5lDlNLyqKOIiKSllRQQurQOpfrzx3AnGVb+NvHn0YdR0Qk7aigHINrzupP5zZ56qWIiMShgnIM2rbI4QfnFfD6ym3MX7M96jgiImlFBeUYXT3iRLq1a8HU2WW4e9RxRETShgrKMWqVl81NFwzk3XU7eH3ltqjjiIikDRWU4/Ct0/vQq2MrpsxeoV6KiEhABeU4tMjJZuLoQt7fsIs5y7ZGHUdEJC2ooBynrw3rRb8urZkyewVVVeqliIiooBynnOwsJo0tYvkne3huyeao44iIRE4FpREuHtKTou5tmTanjIrKqqjjiIhESgWlEbKzjNvGFrGmfB9/XrQp6jgiIpFSQWmkL598AoN7tWf63DIOV6iXIiLNlwpKI5kZk8cNYv2OA/xpwfqo44iIREYFJQFGFeVz2omduHvuKg4eqYw6johIJCIpKGZ2p5ktN7PFZva0mXWso11HM3siaLvMzM4Mlnc2sxIzWxn87JTaPfhcTiaPK+KT3QeZ9c7HUUYREYlMVD2UEmCwuw8ByoCf1NFuOvCiu58EnAIsC5b/GJjr7oXA3GA+UiMLujKyoAv3zVvFvkMVUccREUm5SAqKu8929+rfuvOB3rXbmFl74FzgoWCdw+6+M3h7PDAzeD0TuCy5icOZPK6IbXsPM/PtdVFHERFJuXQ4h3It8EKc5QOAcuARM/ubmf3WzNoE73V3980Awc9udX24mV1nZqVmVlpeXp7o7Ec57cTOnD8onwdeXcPug0eSui0RkXSTtIJiZnPMbGmcaXyNNj8FKoBZcT4iBxgG3OfupwL7OI5DW+4+w92L3b04Pz//OPcmvMnjBrHrwBEeen1t0rclIpJOcpL1we4+pr73zWwCcDEw2uMP2bsB2ODu7wTzT/D3grLFzHq4+2Yz6wGkzQiNg3t14CuDT+ChN9Zyzch+dGqTF3UkEZGUiOoqrwuB24FL3X1/vDbu/gmw3swGBYtGAx8Gr58BJgSvJwB/SWLcYzZpbBH7Dldw/2uro44iIpIyUZ1DuQdoB5SY2SIzux/AzHqa2fM12t0MzDKzxcBQ4BfB8l8CY81sJTA2mE8bRd3bMf6Unsx8ax1b9xyMOo6ISEok7ZBXfdx9YB3LNwEX1ZhfBBTHabedWI8lbU0cU8RfF2/m3ldW86+Xnhx1HBGRpEuHq7yapP5d23D5sN784Z2P2bTzQNRxRESSTgUliW4ZUwjA3S+vijiJiEjyqaAkUa+OrbhyeB/+VLqej7bvizqOiEhSqaAk2Q/PH0h2ljF9zsqoo4iIJJUKSpJ1a9+SCSP78fSijazcsifqOCIiSaOCkgLXnzuA1rnZ/Fq9FBFpwlRQUqBL2xZce3Z/nluymQ827Yo6johIUqigpMg/njOA9i1zmFZSFnUUEZGkUEFJkQ6tcrn+vALmLNvKwo8/jTqOiEjChS4oNYaOl+N0zch+dG6Tx9TZ6qWISNPTYEExs5Fm9iHB0xLN7BQzuzfpyZqgNi1yuHFUAW+s2sbbq7dHHUdEJKHC9FCmAV8GtgO4+/vEnqQox+HqESfSrV0LppasIP6o/SIimSnUIS93X19rUWUSsjQLLXOzufmCgby37lNeW7kt6jgiIgkTpqCsN7ORgJtZnpn9iODwlxyfb53el14dWzFltnopItJ0hCkoNwA/BHoRe4ri0GBejlNeThYTxxSyeMMuZn+4Jeo4IiIJ0WBBcfdt7v5td+/u7t3c/ergeSTSCF87tRf9u7Zh6uwyqqrUSxGRzNfgA7bM7BHgc7/x3P3apCRqJnKys7h1TCETH1vEs0s2c+kpPaOOJCLSKGEOeT0LPBdMc4H2wN5khmouLhnSk0Hd2/HrkjIqKquijiMi0ihhDnk9WWOaBXwTGJz8aE1fVpYxaWwRa7bt4+m/bYw6johIoxzP0CuFQN9EB2muvnxyd77UqwPT567kcIV6KSKSucLcKb/HzHZX/wT+Ctye/GjNg5kxeVwRGz49wOOltW/3ERHJHA2elHf3dqkI0pydV5RP8YmduOfllXzjtN60zM2OOpKIyDGrs4diZsPqm1IZsqmL9VIGsWX3If5n/kdRxxEROS719VCm1POeAxckOEuzdmZBF84a2IX75q3myuF9adOiwc6jiEhaqbOH4u7n1zM1qpiY2Z1mttzMFpvZ02bWsY52Hc3siaDtMjM7M1j+r2a20cwWBdNFjcmTLiaPG8T2fYf53Vvroo4iInLMQl3lZWaDzeybZvbd6qmR2y0BBrv7EKAM+Ekd7aYDL7r7ScApHD2G2DR3HxpMzzcyT1oY1rcTo0/qxgOvrmbXgSNRxxEROSZhrvL6GXB3MJ0P3AFc2piNuvtsd68IZucDveNstz2xYfIfCtY57O47G7PdTDBpbBG7D1bw0Otroo4iInJMwvRQLgdGA5+4+/eI9RRaJDDDtcALcZYPAMqBR8zsb2b221pPjbwpOGT2sJl1SmCeSA3u1YGLvnQCD72xlh37DkcdR0QktDAF5YC7VwEVQa9hK7Ff9vUyszlmtjTONL5Gm58CFcCsOB+RAwwD7nP3U4F9wI+D9+4DCoiNfLyZei4gMLPrzKzUzErLy8tD7G70Jo0pYv+RSh54dXXUUUREQgtzKVFpcNL8QWABsXG83m1oJXcfU9/7ZjYBuBgY7fEfCrIB2ODu7wTzTxAUFHf/bMx3M3uQ2HhjdeWYAcwAKC4uzohhfQu7t+Oyob2Y+fY6vn92f7q1bxl1JBGRBoUZy+tGd9/p7vcDY4EJwaGv42ZmFxK72/5Sd99fx3Y/IfZwr0HBotHAh8H6PWo0/SqwtDF50tHE0YUcqXTunadeiohkhjAn5f9iZleZWRt3X+fuixOw3XuAdkBJcNnv/cG2eppZzSu2bgZmmdliYoe3fhEsv8PMlgTLzwcmJSBTWunXtQ3fLO7NH975mI07D0QdR0SkQdbQI2jN7DzgW8A/EDvU9TjwrLsfTH68xCouLvbS0tKoY4S2cecBzr9zHl8b1otffn1I1HFEpJkyswXuXtxQuzCHvF519xuJnYifQWz4+q2NjygN6dWxFVed0Zc/LdjAum37oo4jIlKvsDc2tgK+Tuz58qcDM5MZSv7uxlEF5GYb0+eujDqKiEi9wpxDeZzYHeoXAL8BCtz95mQHk5hu7Vsy4cx+/HnRRlZu2RN1HBGROoXpoTxCrIjc4O4vB/ekSApdf14BbfJymDanLOooIiJ1CnMO5UV3r0xFGImvc5s8rj27P88v+YSlG3dFHUdEJK7jeQSwROD7Z/enQ6tcppaolyIi6UkFJUN0aJXLdecO4OXlW1nw0adRxxER+ZywV3n1MrORZnZu9ZTsYPJ514zsR5c2eUwtWRF1FBGRz2lwLC8z+y9iNzZ+CFSfS3HgtSTmkjjatMjhB6MK+Plzy3hr9TZGFnSNOpKIyGfC9FAuAwa5+0XufkkwNep5KHL8rh5xIt3bt2Dq7DIaGuVARCSVwhSUNUBusoNIOC1zs7n5gkJKP/qUV8syYzh+EWkewhSU/cAiM3vAzO6qnpIdTOr2zeI+9O7UiinqpYhIGglTUJ4B/h14i9jzUKoniUheThYTRxeyZOMuXvpgS8MriIikQIMn5d19ppnlAUXBohXufiS5saQhXz21F/fNW83UkhWM/WJ3srMs6kgi0syFGctrFLCS2Dhe9wJlumw4ejnZWdw6toiyLXt5dvGmqOOIiIQ65DUFGOfu57n7ucCXgWnJjSVhXPylHpx0Qjt+PWclFZUaYk1EohWmoOS6+2d30rl7GbrqKy1kZRm3jS1i7bZ9PPW3jVHHEZFmLkxBKTWzh8xsVDA9iE7Kp42xX+zOkN4dmD5nJYcqNIaniEQnTEH5AfABcAswkdgd8zckM5SEZ2ZMHjeIjTsP8Mf31kcdR0SasTBXeR0CpgaTpKFzC7tyer9O3P3yKr5R3IeWudlRRxKRZqjOHoqZ/TH4ucTMFteeUhdRGlLdS9m65xD/M/+jqOOISDNVXw9lYvDz4lQEkcYZMaALZw/syr3zVnPl8L60adFg51NEJKHq7KG4++bg5Y3u/lHNCbgxNfHkWEweV8SOfYf53Vvroo4iIs1QmJPyY+Ms+0qig0jjndq3E2O+0I0HXl3NrgMazEBEUqu+cyg/MLMlwKBa50/WAjqHkqYmjS1i98EKfvv6mqijiEgzU18P5Q/AJcQGh7ykxnSau1/dmI2a2Z1mtjwoUE+bWcc4bQaZ2aIa024zuzV4r7OZlZjZyuBnp8bkaUpO7tmBf/hSDx5+Yy3b9x6KOo6INCP1nUPZ5e7r3P3K4LzJAWJPamxrZn0bud0SYLC7DwHKgJ/E2f4Kdx/q7kOB04gNo/908PaPgbnuXgjMDeYlMGlsIQeOVPLAa+qliEjqhBkc8hIzWwmsBV4F1gEvNGaj7j7b3SuC2flA7wZWGQ2sDgobwHhgZvB6JrGnSkpgYLd2XHZqL2a+tY6tuw9GHUdEmokwJ+V/DowAyty9P7Ff7m8mMMO1NFygrgAerTHfvfoqtOBntwTmaRImji6kssr5zSuroo4iIs1EmIJyxN23A1lmluXurwBDG1rJzOaY2dI40/gabX4KVACz6vmcPOBS4E8hssZb/zozKzWz0vLy5vPI3BO7tOEbxX34w7sfs+HT/VHHEZFmIExB2WlmbYHXgFlmNp1YEaiXu49x98Fxpr8AmNkEYjdNftvrf47tV4CF7l7z0YRbzKxH8Dk9gK315Jjh7sXuXpyfn9/gzjYlN18wEMO4e656KSKSfGEKynhiJ8QnAS8Cq2nk3fNmdiFwO3Cpuzf05/OVHH24C2JXnk0IXk8A/tKYPE1Vz46tuOqMvjyxcANrt+2LOo6INHFhCsq/uHuVu1e4+0x3v4tYMWiMe4B2QElwSfD9AGbW08yer25kZq2J3Vj5VK31fwmMDS4WGBvMSxw3nl9AbrYxfU5Z1FFEpImL5E55dx/o7n2qLwt29xuC5Zvc/aIa7fa7exd331Vr/e3uPtrdC4OfOxqTpynr1q4l14zsz1/e30TZlj1RxxGRJizMnfInxblTfknqIkpjXX/uANrk5TCtRL0UEUme+oak/QOxy3n/k6NvHNyjHkFm6dQmj++f3Z/pc1eydOMuBvfqEHUkEWmCGrxTHpgO7Kgx0vARMzsjVQElMb5/Tn86tMplyuwVUUcRkSYqzDmU+4C9Neb3Bcskg7Rvmcv15w3glRXlLPhIHUwRSbwwBcVq3ifi7lWEeHSwpJ9rRvaja9s8pszWuRQRSbwwBWWNmd1iZrnBNBHQqIMZqHVeDjeOGshbq7fz1qptUccRkSYmTEG5ARgJbAQ2AGcA1yUzlCTPVWf05YT2LZlSUkb9AxSIiBybBguKu2919yvcvZu7d3f3q9y9zqFOJL21zM3m5tEDWfDRp8wraz5jm4lI8oUZvr7IzOaa2dJgfoiZ/Z/kR5Nk+cZpfejTuRVTZq9QL0VEEibMIa8HiT0A6wiAuy8mNpy8ZKi8nCwmji5i6cbdvPTBJ1HHEZEmIkxBae3u79Za1uBow5LeLhvakwH5bZhaUkZllXopItJ4YQrKNjMrIPb4X8zscmBzUlNJ0uVkZ3Hb2CLKtuzl2cWboo4jIk1AmILyQ+ABYmN6bQRuJXbll2S4iwb34KQT2jGtpIyKyqqo44hIhgtzldcadx8D5AMnufvZNZ7tLhksK8uYPG4Q67bv56mFG6OOIyIZLsxVXl3M7C7gdWCemU03sy7JjyapMOYL3Tildwemz13JoYrKqOOISAYLc8jrMaAc+DpwefD68WSGktQxi/VSNu48wOPvrY86johksDAFpbO7/7u7rw2mnwMdkx1MUuecwq4M79eZu19exYHD6qWIyPEJU1BeMbMrzCwrmL4JPJfsYJI6sV5KEeV7DvE/83V6TESOT5iCcj2xh20dCqbHgNvMbI+Z7U5mOEmdMwZ04ZzCrtz36mr2HtJtRiJy7MJc5dXO3bPcPTeYsoJl7dy9fSpCSmpMHjeIHfsO87s310YdRUQyUJirvL5faz7bzH6WvEgSlaF9OjLmC9154LU17Np/JOo4IpJhwhzyGm1mz5tZDzP7EjAfaJfkXBKR28YWsedgBQ++rkfeiMixCXPI6ypgJrCE2Mn4W939R8kOJtH4Ys/2/MOQHjz85lq27z0UdRwRySBhDnkVAhOBJ4F1wHfMrHWSc0mEJo0p4uCRSu5/dXXUUUQkg4Q55PVX4F/c/XrgPGAl8F5SU0mkBnZry1dP7c1/v/0RW3YfjDqOiGSIMAVluLvPAfCYKcBljdmomd1pZsvNbLGZPW1mn7tR0swGmdmiGtNuM7s1eO9fzWxjjfcuakwe+byJowuprHJ+88qqqKOISIYIU1AqzOz/mtmD8NkhsEGN3G4JMNjdhwBlxB7gdRR3X+HuQ919KHAasB94ukaTadXvu/vzjcwjtfTt0ppvnt6HR9/9mPU79kcdR0QyQJiC8gixGxrPDOY3AD9vzEbdfba7V989Nx/o3cAqo4HVGuU4tW6+YCBmxt0vr4w6iohkgDAFpcDd7+DvjwA+AFgCM1wLvNBAmyuAR2stuyk4ZPawmXWqa0Uzu87MSs2stLy8vLFZm5UeHVrx7TP68uTCjawp3xt1HBFJc2EKymEza8Xfn9hYQKzHUi8zm2NmS+NM42u0+SmxxwnPqudz8oBLgT/VWHwfUAAMJfb0yCl1re/uM9y92N2L8/PzG4ottdw4aiB52VlMn6teiojULydEm58BLwJ9zGwWcBZwTUMrBQ/lqpOZTQAuBka7e30PNf8KsNDdt9T47M9eB+d2nm0ojxyf/HYtuOasftz/6mpuHDWQQSfonlYRiS/MjY0lwNeIFZFHgWJ3n9eYjZrZhcDtwKXu3tAZ3yupdbjLzHrUmP0qsLQxeaR+1587gLZ5OUwrKYs6ioiksTCHvHD37e7+nLs/6+7bErDde4gN31ISXPZ7P4CZ9TSzz67YCm6gHAs8VWv9O8xsiZktBs4HJiUgk9ShY+s8vn9Of1784BOWbNgVdRwRSVNW/9GmpqW4uNhLS0ujjpGRdh88wrl3vMLQPh353feGRx1HRFLIzBa4e3FD7UL1UETat8zl+nMLmLeinNJ1O6KOIyJpKMxYXr8ys5NTEUbS24SRJ9K1bQumzNa5FBH5vDA9lOXADDN7x8xuMLMOyQ4l6al1Xg4/PL+At9ds561ViTiVJiJNSZirvH7r7mcB3wX6AYvN7A9mdn6yw0n6uXJ4X3p0aMmvZq+gOZ1/E5GGhTqHYmbZwEnBtA14n9hz5R9LYjZJQy1zs7n5gkIWfryTV1ZsjTqOiKSRMOdQpgIrgIuAX7j7ae7+X+5+CXBqsgNK+vlGcW/6dm7NlNllVFWplyIiMWF6KEuBIe5+vbu/W+s9XT/aDOVmZzFxdCEfbNrNSx98EnUcEUkTYQrKIuAkMxtWYyowsxx3111uzdRlp/aiIL8NU0vKqFQvRUQIV1DuJTbE/AzgQeBt4DGgzMzGJTGbpLHsLOO2sYNYuXUvf31/U9RxRCQNhCko64BTgxF7TyN23mQpMAa4I4nZJM19ZfAJfKFHe349p4wjlVVRxxGRiIUpKCe5+wfVM+7+IbECsyZ5sSQTZGUZk8cWsW77fp5csCHqOCISsTAFpczM7jOz84Lp3mBZC4KHbknzNfoL3TilT0fumruSQxWVUccRkQiFKSgTgFXArcRG9V1DbCj7I8RG+pVmzMz40bgiNu06yGPvro86johEqN4HbAU3ND7o7lcT/6mIei6scPbArpzRvzP3vLKKbxb3oVVedtSRRCQC9fZQ3L0SyA8ewysSl5kxedwgyvcc4vfz10UdR0QiEuYRwOuAN83sGWBf9UJ3n5qsUJJ5hvfvzLlF+dw3bzVXnXEibVuE+aclIk1JmHMom4g9sz2L2FMWqyeRo0weW8Sn+4/w8Btro44iIhFo8M9Id/83ADNr4+77GmovzdcpfToy9ovdefC1NXz3zBPp2FpHSkWakzCDQ55pZh8Cy4L5U4JLh0U+57axRew9XMGDr+s2JZHmJswhr18DXwa2A7j7+8C5yQwlmesLPdpz8ZCePPLmOrbtPRR1HBFJoVDPQ3H32jcY6A42qdOtYwo5eKSS++etjjqKiKRQmIKy3sxGAm5meWb2I4LDXyLxFOS35WvDevP7+R/xya6DUccRkRQJU1BuAH4I9AI2AEODeZE6TRxdSGWVc88rK6OOIiIpEuaZ8tvc/dvu3t3du7n71e6+PRXhJHP16dyab53eh8ffW8/6HfujjiMiKRDmKq98M/tnM5thZg9XT43dsJndaWbLzWyxmT1tZh3raDfJzD4ws6Vm9qiZtQyW9zezd8xspZk9rrv5089NFwzEzLhrrnopIs1BmENefwE6AHOA52pMjVUCDHb3IUAZ8JPaDcysF3ALUOzug4Fs4Irg7f8Cprl7IfAp8P0EZJIE6tGhFd8ZcSJPLtzAmnIN+ybS1IUpKK3d/XZ3/6O7P1k9NXbD7j7b3SuC2flA7zqa5gCtzCwHaA1sMjMDLgCeCNrMBC5rbCZJvB+MKqBFTja/nqNeikhTF6agPGtmFyU5x7XAC7UXuvtG4FfAx8BmYJe7zwa6ADtrFKQNxC4akDTTtW0LvndWP/66eBPLP9kddRwRSaIwBWUisaJy0Mx2m9keMwv1m8HM5gTnPmpP42u0+SlQAcyKs34nYDzQH+gJtDGzqwGLszmvI8N1ZlZqZqXl5eVhYkuCXXfuANrm5TB1dlnUUUQkicKM5XXcA0G6+5j63jezCcDFwGh3j1cQxgBr3b08aP8UMJJY8eloZjlBL6U3sUEs42WYAcwAKC4ujlt0JLk6ts7jH88ZwLQ5ZSzesJMhveNefyEiGS7MVV5mZleb2f8N5vuY2fDGbtjMLgRuBy5197quK/0YGGFmrYPzJqOBZUHxeQW4PGg3gdjFA5Kmrj27H51a5zJFvRSRJivMIa97gTOBq4L5vcBvErDte4gNg19iZovM7H4AM+tpZs8DuPs7xE68LwSWBHlnBOvfDtxmZquInVN5KAGZJEnatczlhvMKeLWsnPfW7Yg6jogkgcU/0lSjgdlCdx9mZn9z91ODZe+7+ykpSZhAxcXFXlpaGnWMZuvA4UrOueMVCvLb8Nh1I4h1OkUk3ZnZAncvbqhdmB7KkeDZ8h58cD5Q1ch80gy1ysvmpvMLeGftDt5arcEWRJqaMAXlLuBpoJuZ/QfwBvCLpKaSJuvKM/rSs0NL7nxpBQ31jkUks4QZy2sW8E/AfxK7F+Qyd/9TsoNJ09QiJ5ubRxeyaP1OXl6+Neo4IpJAYZ+Hstzdf+Pu97i7hq6XRrn8tN6c2KU1U2aXUVWlXopIUxGqoIgkUm52FreOKeTDzbt58YNPoo4jIgmigiKRuPSUXgzs1papJWVUqpci0iSooEgksrOM28YWsWrrXp55f2PUcUQkAVRQJDIXnnwCX+zRnmklKzlSqSvRRTKdCopEJivLmDyuiI937OeJBRuijiMijaSCIpG64KRunNq3I3fNXcnBI5VRxxGRRlBBkUiZGT8aN4jNuw7y2LsfRx1HRBpBBUUiN7KgCyMGdOaeV1Zz4LB6KSKZSgVFImdmTB43iG17D/Hfb6+LOo6IHCcVFEkLp/frzHlF+dz/6mr2HDwSdRwROQ4qKJI2Jo8r4tP9R3j4jXVRRxGR46CCImljSO+OfPnk7vz29TXs3H846jgicoxUUCStTBpbxN7DFcx4bU3UUUTkGKmgSFo56YT2XDKkJ4+8uY5tew9FHUdEjoEKiqSdW8cUcqiikvvmrY46iogcAxUUSTsD8tvy9WG9+f38j/hk18Go44hISCookpZuGV2Iu3P3yyujjiIiIamgSFrq07k1V5zel8ffW8/6HfujjiMiIaigSNq66YKBZGcZ0+eqlyKSCVRQJG11b9+S74w4kacWbmB1+d6o44hIA1RQJK3dMKqAlrnZ/HqOeiki6S6SgmJmd5rZcjNbbGZPm1nHOtpNMrMPzGypmT1qZi2D5b8zs7VmtiiYhqZ2DyRVurZtwffO6sdf39/Ess27o44jIvWIqodSAgx29yFAGfCT2g3MrBdwC1Ds7oOBbOCKGk3+t7sPDaZFqQgt0bjunALatcxhaklZ1FFEpB6RFBR3n+3uFcHsfKB3HU1zgFZmlgO0BjalIp+klw6tc7nunAGUfLiF99fvjDqOiNQhHc6hXAu8UHuhu28EfgV8DGwGdrn77BpN/iM4ZDbNzFqkJqpE5Xtn96dT61ymqJcikraSVlDMbE5w7qP2NL5Gm58CFcCsOOt3AsYD/YGeQBszuzp4+yfAScDpQGfg9npyXGdmpWZWWl5enrD9k9Rq2yKHH4wq4LWyct5duyPqOCISR9IKiruPcffBcaa/AJjZBOBi4Nvu7nE+Ygyw1t3L3f0I8BQwMvjszR5zCHgFjd3QAAAIGUlEQVQEGF5PjhnuXuzuxfn5+YneTUmh74zoR367Fvxq9gri/5MRkShFdZXXhcR6FZe6e123QX8MjDCz1mZmwGhgWbB+j+CnAZcBS5OfWqLWKi+bm84fyLtrd/Dmqu1RxxGRWqI6h3IP0A4oCS77vR/AzHqa2fMA7v4O8ASwEFgSZJ0RrD/LzJYEy7sCP09xfonIFcP70KtjK+5UL0Uk7eREsVF3H1jH8k3ARTXmfwb8LE67C5KXTtJZi5xsbhk9kNufXMLcZVsZ88XuUUcSkUA6XOUlcky+Nqw3/bq0ZkpJGVVV6qWIpItIeigijZGbncWtY4q49fFFDPzp82RnGVlmZGcZ2WZkZVmNZZBlNd7PMrKMo9ep8TrLqGN57LNqL6treXatDNlmWJzl8bYVa0sdGWqtd9TnUk+Guj839t+MOMss6q9aMowKimSkS07pyY59h9m+7xBVDlVVTmWVU+kee+1OZRWfvf7sZ51tnarqn1VwqKKKqgY+K9Y29pnxlseW8dmyTPS54hunYFcvq25jBk2pFMWu/cl8v/jqlxjev3NSt6GCIhkpO8u49uz+Ucc4JlXVRatWoYlfqKjR9ujlny+AtYpancWydobqtsQvhl6riNZXsD9bHsvdZDShXWnTIjvp21BBEUmRrCwjC9P/dNJk6aS8iIgkhAqKiIgkhAqKiIgkhAqKiIgkhAqKiIgkhAqKiIgkhAqKiIgkhAqKiIgkhDWnIcDNrBz46DhX7wpsS2CcKGlf0k9T2Q/QvqSrxuzLie7e4BMKm1VBaQwzK3X34qhzJIL2Jf00lf0A7Uu6SsW+6JCXiIgkhAqKiIgkhApKeDMabpIxtC/pp6nsB2hf0lXS90XnUEREJCHUQxERkYRQQanFzC40sxVmtsrMfhzn/RZm9njw/jtm1i/1KcMJsS/XmFm5mS0Kpn+MImdDzOxhM9tqZkvreN/M7K5gPxeb2bBUZwwjxH6MMrNdNb6Pf0l1xrDMrI+ZvWJmy8zsAzObGKdNpnwvYfYl7b8bM2tpZu+a2fvBfvxbnDbJ/f3l7pqCCcgGVgMDgDzgfeCLtdrcCNwfvL4CeDzq3I3Yl2uAe6LOGmJfzgWGAUvreP8i4AViT54dAbwTdebj3I9RwLNR5wy5Lz2AYcHrdkBZnH9fmfK9hNmXtP9ugv/ObYPXucA7wIhabZL6+0s9lKMNB1a5+xp3Pww8Boyv1WY8MDN4/QQw2tLzodNh9iUjuPtrwI56mowH/ttj5gMdzaxHatKFF2I/Moa7b3b3hcHrPcAyoFetZpnyvYTZl7QX/HfeG8zmBlPtk+RJ/f2lgnK0XsD6GvMb+Pw/rM/auHsFsAvokpJ0xybMvgB8PTgc8YSZ9UlNtIQLu6+Z4MzgkMULZnZy1GHCCA6bnErsL+KaMu57qWdfIAO+GzPLNrNFwFagxN3r/E6S8ftLBeVo8Sp17Qofpk06CJPzr0A/dx8CzOHvf7lkmkz5ThqykNgQF6cAdwN/jjhPg8ysLfAkcKu77679dpxV0vZ7aWBfMuK7cfdKdx8K9AaGm9ngWk2S+p2ooBxtA1Dzr/TewKa62phZDtCB9DyM0eC+uPt2dz8UzD4InJaibIkW5ntLe+6+u/qQhbs/D+SaWdeIY9XJzHKJ/QKe5e5PxWmSMd9LQ/uSad+Nu+8E5gEX1norqb+/VFCO9h5QaGb9zSyP2EmrZ2q1eQaYELy+HHjZgzNcaabBfal1PPtSYseOM9EzwHeDq4pGALvcfXPUoY6VmZ1QfTzbzIYT+/9ze7Sp4gtyPgQsc/epdTTLiO8lzL5kwndjZvlm1jF43QoYAyyv1Sypv79yEvVBTYG7V5jZTcBLxK6SetjdPzCz/weUuvszxP7h/d7MVhGr7FdEl7huIfflFjO7FKggti/XRBa4Hmb2KLGrbLqa2QbgZ8ROOOLu9wPPE7uiaBWwH/heNEnrF2I/Lgd+YGYVwAHgijT9YwXgLOA7wJLgmD3APwN9IbO+F8LtSyZ8Nz2AmWaWTazg/dHdn03l7y/dKS8iIgmhQ14iIpIQKigiIpIQKigiIpIQKigiIpIQKigiIpIQKigix8HM3gp+9jOzqxL82f8cb1si6U6XDYs0gpmNAn7k7hcfwzrZ7l5Zz/t73b1tIvKJpJJ6KCLHwcyqR3X9JXBO8IyMScHgfHea2XvBoJvXB+1HBc/c+AOwJFj2ZzNbEDy74rpg2S+BVsHnzaq5reCO8zvNbKmZLTGzb9X47HnBAJ/LzWxWmo6ALU2c7pQXaZwfU6OHEhSGXe5+upm1AN40s9lB2+HAYHdfG8xf6+47gmEy3jOzJ939x2Z2UzDAX21fA4YCpwBdg3VeC947FTiZ2FhZbxK7+/uNxO+uSN3UQxFJrHHExq9aRGwI9C5AYfDeuzWKCcSGvnkfmE9swL5C6nc28GgwouwW4FXg9BqfvcHdq4BFQL+E7I3IMVAPRSSxDLjZ3V86amHsXMu+WvNjgDPdfb+ZzQNahvjsuhyq8boS/b8tEVAPRaRx9hB7bGy1l4gNIpgLYGZFZtYmznodgE+DYnISsUfkVjtSvX4trwHfCs7T5BN7pPC7CdkLkQTQXzEijbMYqAgOXf0OmE7scNPC4MR4OXBZnPVeBG4ws8XACmKHvarNABab2UJ3/3aN5U8DZwLvE3so0j+5+ydBQRKJnC4bFhGRhNAhLxERSQgVFBERSQgVFBERSQgVFBERSQgVFBERSQgVFBERSQgVFBERSQgVFBERSYj/DxUmhi71gb99AAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "import scipy.optimize\n",
    "import matplotlib.pyplot as plt\n",
    "\n",
    "def cost(phi):\n",
    "    return ExpectVal(M, TwoQubitPQC(phi))\n",
    "\n",
    "cost_val = [] #コスト関数の変化を保存するための関数\n",
    "\n",
    "#この関数がiteration ごとに呼ばれる。\n",
    "def callback(phi):\n",
    "    global cost_val\n",
    "    cost_val.append(cost(phi))\n",
    "\n",
    "init = np.random.rand(n_param)\n",
    "callback(init)\n",
    "res = scipy.optimize.minimize(cost, init, \n",
    "                              method='Powell',\n",
    "                              callback=callback)\n",
    "plt.plot(cost_val)\n",
    "plt.xlabel(\"iteration\")\n",
    "plt.ylabel(\"energy expectation value\")\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ハミルトニアンを対角化して得られた厳密なエネルギーと比べることで、VQE によって算出された値が正しいか検証してみよう。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2.8626207640766816\n",
      "-2.8623984117519257\n"
     ]
    }
   ],
   "source": [
    "import scipy.linalg\n",
    "l, P = scipy.linalg.eigh(M)\n",
    "print(l[0]) #最小固有値\n",
    "print(cost(res.x)) #VQEの結果"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Powell法で算出した固有値と一致はしていないものの、小数第3位まで同じであるので、殆ど正しいと言っていいだろう。　\n",
    "\n",
    "次に回路の出力にノイズが存在するケースでも検証してみよう。NISQでは出力にエラー（ノイズ）がのることが避けられないため、ノイズありでもアルゴリズムが動くのか・どの程度のノイズまでなら耐えられるのかといった検証は非常に重要である。"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "def cost(phi):\n",
    "    return ExpectVal(M,TwoQubitPQC(phi))+np.random.normal(0,0.01)\n",
    "\n",
    "def callback(phi):\n",
    "    global cost_val\n",
    "    cost_val.append(cost(phi))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY4AAAEKCAYAAAAFJbKyAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADl0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uIDMuMC4yLCBodHRwOi8vbWF0cGxvdGxpYi5vcmcvOIA7rQAAIABJREFUeJzt3Xl0XPV99/H3V7ssy5IXeZGlsdkdMMY2QgJCCAlOIMQYEkgwMi0tfUpo0oa0JycJzdPk9Hm606RNmiYU0uZJWwvMEsK+mBBCSINt2RjbYGPAgBd5kW0sb5JtSd/nj7nCsjyS7kiauSPN53XOPZ7lN3M/Ghh9dX/fu5i7IyIiElZO1AFERGR4UeEQEZGkqHCIiEhSVDhERCQpKhwiIpIUFQ4REUmKCoeIiCRFhUNERJKiwiEiIknJizpAKkyYMMGnT58edQwRkWFj5cqVu929IszYEVk4pk+fTmNjY9QxRESGDTN7L+xYTVWJiEhSVDhERCQpKhwiIpIUFQ4REUmKCoeIiCRFhUNERJKiwiEiIklR4Qi0Hevgnhc38du390QdRUQko6lwBHJzjLt/vYl7fr0p6igiIhlNhSOQn5vDwguqeeGNXWzb1xp1HBGRjKXC0c0NF1TjwJLlm6OOIiKSsVQ4uqkaO4rLzqzgvhVbONbRGXUcEZGMpMLRQ33dNHYdOMIv1u+KOoqISEZS4ejhY2dVMHlMEQ2arhIRSUiFo4e83BwW1lbz6zeb2bL3cNRxREQyjgpHAjdcUI0B92qrQ0TkJCocCUwpK+bjMyZxf+MWjrarSS4i0p0KRy8W1cXYffAoS1/fGXUUEZGMosLRi0vPrGBqeTENy0NfTVFEJCuocPQiN8e4sbaa37y1h3d3H4o6johIxlDh6MPna6rJyzE1yUVEulHh6MPEMUXM+9AkHli5lSPtHVHHERHJCCoc/aivi7H30FGeXrcj6igiIhlBhaMfl5w+gdi4UTQs03SViAhEVDjM7E4z22Bma8zsYTMr72Nsrpm9YmaPpzNjl5wc48baGMve2ctbuw5GEUFEJKNEtcWxFJjp7rOAjcAdfYy9HVifllS9+FxNFfm5pq0OEREiKhzu/qy7twd3XwaqEo0zsyrg08CP05UtkQmjC/nkOZN5aNVW2o6pSS4i2S0Tehy3AE/18tw/A18D+j3vh5ndamaNZtbY3Nw8lPkAWFQbo6X1GE+u3T7k7y0iMpykrHCY2XNmti7Bck23Md8E2oHFCV4/H9jl7ivDrM/d73b3GnevqaioGLKfo8tFp43n1Aklmq4SkayXl6o3dvd5fT1vZjcD84HL3d0TDPkwsMDMrgKKgDFm9t/uftPQp+2fWbxJ/tdPrueNHQc4a3JpFDFERCIX1V5VVwJfBxa4e8KLXrj7He5e5e7TgYXA81EVjS7XnV9FQW4ODct0/ioRyV5R9Th+AJQCS81stZndBWBmlWb2ZESZ+jWupIBPnTuZn72yjdajapKLSHaKaq+q09292t1nB8ttweNN7n5VgvEvuPv89Cc9WX1tjANt7Ty2pinqKCIikciEvaqGldpTxnH6xNFqkotI1lLhSJKZUV8bY/WWfbzW1BJ1HBGRtFPhGIDr5lZRmJejrQ4RyUoqHANQNiqfT8+awiOrmzh0pL3/F4iIjCAqHAO0qC7GwSPtPPqqmuQikl1UOAZobmwsMyaXarpKRLKOCscAmRn1dTHWbmthzdZ9UccREUkbFY5BuHbOVIrzc7XVISJZRYVjEMYU5XP1eVN49NUmDrQdizqOiEhaqHAM0qK6aRw+2sHPV6tJLiLZQYVjkGZVlXFO5Rgalm0m8Ul+RURGltCFw8xKUhlkuOpqkq/fvp9XtqhJLiIjX7+Fw8wuNrPXCa77bWbnmdkPU55sGLlm9lRKCtQkF5HsEGaL45+AK4A9AO7+KnBpKkMNN6ML81gweyqPr2mipVVNchEZ2UJNVbn7lh4P6WIUPSyqi9F2rJOHV22NOoqISEqFKRxbzOxiwM2swMy+SjBtJcfNnFrGeVVlNCxXk1xERrYwheM24EvAVGArMDu4Lz3U18XYuPMgje+9H3UUEZGU6bdwuPtud1/k7pPcfaK73+Tue9IRbri5+rxKSgvz1CQXkREtr78BZvYT4KS5F3e/JSWJhrFRBXlcO2cqSxq38K35ZzO2pCDqSCIiQy7MVNXjwBPB8gtgDHAwlaGGs/q6GEfbO3lITXIRGaH63eJw94e63zeze4HnUpZomPvQlDHMjZXTsHwzf3DJKZhZ1JFERIbUQE45cgYQG+ogI0l93TQ2NR/i5U17o44iIjLkwhw5fsDM9nf9CzwGfD310Yav+bOmMKYoj4blapKLyMgTZqqqNB1BRpKi/Fw+O7eKxcveY8/Bsxk/ujDqSCIiQ6bXLQ4zm9vXks6Qw9GiuhjHOpwHV6pJLiIjS19bHN/p4zkHPj7EWUaUMyaVUjt9HA3LN/OHHzmVnBw1yUVkZOi1cLj7x9IZZCSqr4vxlSWr+Z+393DJGROijiMiMiT67XEAmNlM4GygqOsxd//PVIUaKa6cOZmxj+XTsPw9FQ4RGTHC7FX1beBfguVjwD8AC1Kca0Qoys/lurlVPPvaTnYdaIs6jojIkAhzHMf1wOXADnf/feA8YFC7CZnZnWa2wczWmNnDZlbey7h3zWytma02s8bBrDMqN9bFaO90HmhUk1xERoYwhaPV3TuBdjMbA+wCTh3kepcCM919FrARuKOPsR9z99nuXjPIdUbitIrRXHTqeO5dvpnOTp1uXUSGvzCFozHYIrgHWAmsApYPZqXu/qy7twd3XwaqBvN+ma6+LsbW91t58c3mqKOIiAxamNOqf9Hd97n7XcAngJuDKauhcgvwVG+rB541s5VmdusQrjOtrjhnMuNLCnS6dREZEcKcVv0RYAnwiLu/G/aNzew5YHKCp77p7o8EY74JtAOLe3mbD7t7k5lNBJaa2QZ3f7GX9d0K3AoQi2XWqbQK8nK4vqaKH//6HXa0tDG5rKj/F4mIZKgwU1XfBS4BXjezB8zsejPr9zefu89z95kJlq6icTMwH1jkvVxr1d2bgn93AQ8DtX2s7253r3H3moqKihA/VnrV18bo6HSWrOh5+XYRkeElzFTVr9z9i8Qb4ncDnyfeIB8wM7uS+IkSF7j74V7GlJhZaddt4JPAusGsN0rTxpfwkTMmsGTFZjrUJBeRYSzUadXNrBi4jvj1xy8AfjrI9f4AKCU+/bTazO4K1lNpZk8GYyYBL5nZq8Sb8U+4+9ODXG+k6mtjNLW08cIbg6q7IiKRCtPjWALUAU8D/wq8EOyeO2DufnovjzcBVwW3NxE/ZmTEmHf2JCpKC2lYtpnLPzQp6jgiIgMSZovjJ8Bp7n6buz8/2KKRzfJzc7ihpppfvrGLbftao44jIjIgYXocT7t7RzrCZIOFtdU4sEQXeRKRYWogl46VQagaO4qPnlnBksYttHdo401Ehh8VjgjU18bYuf8Iv9igJrmIDD9hT6s+FZjWfXxvB+JJ/z4+YyKTxxTRsGwzV5yT6BhJEZHMFWavqr8HbgBeB7p6HQ6ocAxQXm4ON1xQzfeff5Mtew9TPW5U1JFEREILM1V1LXCWu1/l7lcHi67HMUgLa6sx4F41yUVkmAlTODYB+akOkm2mlBXz8RkTub9xK8fUJBeRYSRM4TgMrDazfzOz73ctqQ6WDerrYuw+eISlr++MOoqISGhhmuOPBosMsY+eOZGp5cU0LNvMVedOiTqOiEgo/RYOd/+pmRUAZwYPveHux1IbKzvk5hgLL6jmO0s38u7uQ0yfUBJ1JBGRfvU7VWVmlwFvEj9P1Q+BjWZ2aYpzZY3PX1BNbo6pSS4iw0aYHsd3gE+6+0fd/VLgCuCfUhsre0waU8S8D03kgZVbOdKuM7uISOYLUzjy3f2NrjvuvhHtZTWk6uumsffQUZ55TU1yEcl8YQpHo5n9u5ldFiz3ACtTHSybfOT0CVSPK6Zh2XtRRxER6VeYwvFHwGvAl4HbiR9BflsqQ2WbnBzjxtoYL2/ay1u7DkYdR0SkT2FOq37E3b/r7p9198+4+z+5+5F0hMsmnzu/mjw1yUVkGOi1cJjZ/cG/a81sTc8lfRGzQ0VpIVecM5mHVm2l7Zia5CKSufo6juP24N/56Qgi8SPJn1i7nafWbeczc6qijiMiklCvWxzuvj24+UV3f6/7AnwxPfGyy0Wnjmf6+FE0LNN0lYhkrjDN8U8keOxTQx1E4k3y+roYK959n407D0QdR0Qkob56HH9kZmuBs3r0N94B1ONIkevPr6YgN0dbHSKSsfra4mgAriZ+gsOruy3nu/tNaciWlcaVFHDlzHiTvPWomuQiknn66nG0uPu77n5j0NdoJX7lv9FmFktbwixUXxfjQFs7j69pijqKiMhJwpzk8GozexN4B/gV8C7wVIpzZbW6U8ZxWkUJDTqmQ0QyUJjm+F8BFwIb3f0U4HLgNylNleXMjPq6abyyeR+vN+2POo6IyAnCFI5j7r4HyDGzHHf/JTA7xbmy3nVzp1KQl0PDcp2/SkQyS5jCsc/MRgMvAovN7HtAe2pjSfmoAuafO4Wfv9LEoSP6uEUkc4QpHNcQv+74nwJPA2+jo8nTor4uxsEj7Tz2qprkIpI5whSOb7l7p7u3u/tP3f37wNdTHUzg/GljOWtSKYt1TIeIZJBIjhw3szvNbENwQOHDZlbey7hyM3swGLvezC4azHqHm3iTPMbabS2s3doSdRwRESDckeMzEhw5vnaQ610KzHT3WcBG4I5exn0PeNrdZwDnAesHud5h59o5UynKV5NcRDJHmCPHH+HkI8cXDWal7v6su3d1fF8GTjoVrJmNAS4F/j14zVF33zeY9Q5HZcX5XD2rkkdWN3Gg7VjUcURE+j9ynPhf/Xu7nRn3mJnVDWGGW0h8QOGpQDPwEzN7xcx+bGYlQ7jeYaO+Lsbhox08slpNchGJXpgex4+A7tczPRQ81icze87M1iVYruk25pvEd+1dnOAt8oC5wI/cfU6w3m/0sb5bzazRzBqbm5tD/FjDx+zqcs6eMobFyzbj7lHHEZEsF6ZwmHf7beXunfR9AaiucfPcfWaC5REAM7uZ+G69izzxb8OtwFZ3Xxbcf5B4IeltfXe7e42711RUVIT4sYaPrib5+u37Wb0l62brRCTDhCkcm8zsy2aWHyy3A5sGs1Izu5L4Lr0L3P1wojHuvgPYYmZnBQ9dDrw+mPUOZ9fMrmRUQa5Oty4ikQtTOG4DLga2Ed8KqANuHeR6fwCUAkvNbLWZ3QVgZpVm9mS3cX9C/Gj1NcRPc/I3g1zvsFValM81syt5bE0TLa1qkotIdMJMOe0CFg7lSt399F4ebwKu6nZ/NVAzlOsezuprp3Hv8i38/JVt3Hzx9KjjiEiWCnNa9TPN7Bdmti64P8vM/nfqo0lP51aVMauqjMXL3lOTXEQiE2aq6h7iB+gdA3D3NQzxFoiEV18bY+POg6x87/2oo4hIlgpTOEa5+/Iej+l0rRG5+rxKRhfmqUkuIpEJUzh2m9lpxC8bi5ldD2xPaSrpVUlhHtfOqeTxtdvZd/ho1HFEJAuFKRxfAv6N+DmrtgFfIb6nlUSkvnYaR9s7eXDl1qijiEgW6rdwuPsmd58HVAAz3P2S4NQjEpGzK8cwJ1ZOw3IdSS4i6Rdmr6rxZvZ94NfAC2b2PTMbn/po0pf62hibmg+x7J29UUcRkSwTZqrqPuInG7wOuD64vSSVoaR/82dVUlqkJrmIpF+YwjHO3f+vu78TLH8FJLzwkqRPcUEu182t4ul1O9hz8EjUcUQki4QpHL80s4VmlhMsnweeSHUw6V99XYyjHWqSi0h6hSkcXyB+UacjwXIf8GdmdsDM9qcynPTtzEmlXDB9LPcu30xnp5rkIpIeYfaqKnX3HHfPD5ac4LFSdx+TjpDSu/q6GO/uOcxvN+2JOoqIZIkwe1X9QY/7uWb27dRFkmR8auYUykflq0kuImkTZqrqcjN70symmNm5xK8RXpriXBJSUX68Sf7MaztoPqAmuYikXpipqnrgp8Ba4k3xr7j7V1MdTMKrr4vR3unc37gl6igikgXCTFWdAdwOPAS8C/yOmY1KcS5JwmkVo7nw1HHct0JNchFJvTBTVY8B33L3LwAfBd4EVqQ0lSStvm4aW/a28uu3dkcdRURGuDCFo9bdnwPwuO8A16Y2liTrinMmMa6kgIZlOo2YiKRWmMLRbmZ/YWb3wAdTV2elNpYkqzAvl8+dX8Vz63exc39b1HFEZAQLUzh+QvzAv4uC+1uBv0pZIhmwG2tjdHQ6S1aoSS4iqROmcJzm7v/A8UvHtgKW0lQyINMnlHDJ6RO4b/lmOtQkF5EUCVM4jppZMcevAHga8S0QyUD1dTGaWtr41cZdUUcRkREqTOH4NvA0UG1mi4FfAF9LaSoZsE+cPYkJowt1JLmIpExefwPcfamZrQIuJD5Fdbu7a5/PDJWfm8Pna6q461dv07Svlcry4qgjicgIE2aLA3ff4+5PuPvjKhqZ78baGA7cpya5iKRAqMIhw0v1uFFcekYFS1Zspr2jM+o4IjLCqHCMUPV1MXbuP8LzG9QkF5GhFeZcVf9oZuekI4wMnctnTGTSmEIalqtJLiJDK8wWxwbgbjNbZma3mVlZqkPJ4OXl5nBDTTW/2tjMlr2Ho44jIiNImNOq/9jdPwz8LjAdWGNmDWb2sVSHk8G5oTaGAfet0FaHiAydUD0OM8sFZgTLbuBV4tcdv28gKzWzO81sg5mtMbOHzaw8wZizzGx1t2W/mX1lIOvLVlPLi/nYWRO5v3Erx9QkF5EhEqbH8V3gDeAq4G/c/Xx3/3t3vxqYM8D1LgVmuvssYCNwR88B7v6Gu89299nA+cBh4OEBri9r1dfFaD5whOde3xl1FBEZIcJscawDZrn7F9x9eY/nageyUnd/1t3bg7svA1X9vORy4G131znDk3TZWROpLCtisY4kF5EhEqZwrAZmmNncbstpZpbn7i1DkOEW4Kl+xiwE7u1rgJndamaNZtbY3Nw8BLFGhtwc44YLYrz01m7e3X0o6jgiMgKEKRw/JL5VcDdwD/Bb4D5go5l9srcXmdlzZrYuwXJNtzHfBNqBxX28TwGwAHigr5Dufre717h7TUVFRYgfK3vccEE1uTnGvWqSi8gQCFM43gXmBL+Uzyfe11gHzAP+obcXufs8d5+ZYHkEwMxuBuYDi9y9r3OAfwpY5e6apB+gyWVFXD5jIg82buVou5rkIjI4YQrHDHd/reuOu79OvJBsGuhKzexK4OvAAnfv7yCDG+lnmkr6V18XY8+hozzz2o6oo4jIMBemcGw0sx+Z2UeD5YfBY4UEF3cagB8ApcDSYFfbuwDMrNLMnuwaZGajgE8APxvgeiRw6RkVVI0tZrGuSS4ig9TvadWBm4EvAl8hflr1l4CvEi8aAzoI0N1P7+XxJuK7/XbdPwyMH8g65EQ5OcaNtTHufOYN3m4+yGkVo6OOJCLDVJ9bHMGBf/e4+3fc/TPufq27/6O7H3b3Tnc/mKacMgQ+V1NFXo5xr3bNFZFB6LNwuHsHUBHs2STD3MTSIj55ziQeXLWVtmMdUccRkWEqzFTVu8BvzOxR4IMDAdz9u6kKJalTXzuNJ9fu4Ol1O7h2ztSo44jIMBSmOd4EPB6MLe22yDB08WnjmTZ+lJrkIjJgYa45/pcAZlbi7jr0eJjLyTHqa2P87VMb2LjzAGdO0t8AIpKcMCc5vMjMXgfWB/fPC3bJlWHq+vOrKMjNoUFNchEZgDBTVf8MXAHsAXD3V4FLUxlKUmv86EKumDmZn63aSutRNclFJDmhrsfh7lt6PKTfNsNcfW2M/W3tPL6mKeooIjLMhCkcW8zsYsDNrMDMvkowbSXD14WnjuPUihJdk1xEkhamcNwGfAmYCmwFZgf3ZRgzizfJX9m8j/Xb90cdR0SGkTDXHN/t7ovcfZK7T3T3m9x9TzrCSWpdf34VBXlqkotIcvrdHdfMKoA/BKZ3H+/ut6QulqRD+agCPn3uFB5+ZRvf+NQMSgrDHA8qItkuzFTVI0AZ8BzwRLdFRoD6uhgHj7Tz2KtqkotIOGH+xBzl7l9PeRKJRM20sZwxcTQNyzezsDYWdRwRGQbCbHE8bmZX9T9MhiMzY1FdjDVbW1i3bSguIS8iI12YwnE78eLRZmb7zeyAmWk3nBHkM3OrKMrPYbGa5CISQpi9qkrdPcfdi9x9THB/TDrCSXqUFeczf1Ylj67exsEj7VHHEZEMF+ZcVWZmN5nZXwT3q82sNvXRJJ3q62IcOtrBz1/ZFnUUEclwYaaqfghcBNQH9w8C/5qyRBKJOdXlzJhcSsOyzbh71HFEJIOFKRx17v4loA3A3d8HdEXAEcbMWHThNF7fvp9Xt6pJLiK9C1M4jgXXHnf44IDAzpSmkkhcO7uSUQW5NOgiTyLShzCF4/vAw8BEM/tr4CXgb1KaSiJRWpTPgvMqefTVJlpaj0UdR0QyVJi9qhYDXwP+FtgOXOvuD6Q6mESjvi5G27FONclFpFehTk7k7huADSnOIhlgVlU5504to2HZZn73ommYWdSRRCTDhLqQk2SX+roYb+w8wKrN70cdRUQykAqHnGTBeZWMLszTkeQikpAKh5ykpDCPa2ZX8via7ew7fDTqOCKSYVQ4JKH6uhhH2zt5aJWa5CJyIhUOSeicyjJmV5fTsOw9HUkuIieIpHCY2Z1mtsHM1pjZw2ZW3su4PzWz18xsnZnda2ZF6c6azerrYrzdfIjl7+yNOoqIZJCotjiWAjPdfRawEbij5wAzmwp8Gahx95lALrAwrSmz3NWzKiktyqNhuZrkInJcJIXD3Z91967zd78MVPUyNA8oNrM8YBSg65umUXFBLp+dM5Wn1u5g7yE1yUUkLhN6HLcAT/V80N23Af8IbCZ+xHqLuz+b5mxZr75uGkc7Onlw5Zaoo4hIhkhZ4TCz54LeRM/lmm5jvgm0A4sTvH4scA1wClAJlJjZTX2s71YzazSzxubm5qH/gbLUWZNLqZk2lnuXb1GTXESAFBYOd5/n7jMTLI8AmNnNwHxgkSf+jTQPeMfdm939GPAz4OI+1ne3u9e4e01FRUUqfqSsVV8X453dh/jt23uijiIiGSCqvaquBL4OLHD3w70M2wxcaGajLH7CpMuB9enKKMddde4UyorzdSS5iADR9Th+AJQCS81stZndBWBmlWb2JIC7LwMeBFYBa4Osd0eUN6sV5edy3dwqnnltB80HjkQdR0QiFtVeVae7e7W7zw6W24LHm9z9qm7jvu3uM4Iprt9xd/3Wikh9XTXtnc4DapKLZL1M2KtKhoHTJ5ZSd8o47lu+hc5ONclFspkKh4RWXxdj897DvPTW7qijiEiEVDgktCtnTmZcSQGLdU1ykaymwiGhFeblcv35VTy3fhc797dFHUdEIqLCIUm5sTZGR6dz/wo1yUWylQqHJOWUCSV8+PTx3LdiCx1qkotkJRUOSVp97TS27WvlxY06tYtINlLhkKR94uxJTBitJrlItlLhkKQV5OXwuZpqnt+wi6Z9rVHHEZE0y4s6gAxPN14Q40cvvM0X/msl51aVMbW8mCllRVSWF1NZVszksiIK8vR3ichIpMIhAxIbP4qvXXkWT6zZzlNrt/P+4WMnPG8GE0YXUhkUkyllxVSWB4WlvJjKsiImjC4kJ8ci+glEZKBsJF5joaamxhsbG6OOkVVaj3bQ1NLK9n1tNO1rpamllaZ9rWxvaWPbvvjjrcc6TnhNfq4xuayIyrLioLh0FZbjxWZMUR7xkyOLSCqZ2Up3rwkzVlscMiSKC3I5rWI0p1WMTvi8u7Pv8LHjxaWllaagyGxvaWX5O3vZsb/tpF18RxfmnVhQyoqZEmyxVJbHp8SK8nPT8SOKSECFQ9LCzBhbUsDYkgLOqSxLOKaj09l1oI2mfW1sD7ZYjheXNtZta2FPgmufTxhd8MFU2JSy4ni/pfx4v6WitJBcTYmJDBkVDskYuTnGlLL4FBWMTTim7VgH21va2L6vlaaW41ss2/a1san5EC+9uZtDR0+cEsvLMSaNKTphCmxqUGS6tmTKivM1JSYSkgqHDCtF+bmcMqGEUyaUJHze3dnf1n5CQdnerdeyavP77GjZzrGOE6fERhXknrBXWPctlq6CoykxkTgVDhlRzIyy4nzKivP50JQxCcd0djq7Dx6JN+2DrZbu/ZYNOw4kvNLhuJKCbsUl2HrpdntiaSF5udoFWUY+FQ7JOjk5xsQxRUwcU8ScXsYcae9gZ0tXcQn6LUGR2bznMC9v2sOBtvYTXpObY0wqLYwXk+7Fpez4bsilRXkYfDAt1jU5ZoamymTYUOEQSaAwL5fY+FHExo/qdcyBtmMn7G7cfTfkNVv38cy6No52dA44Q1cd6a3QxO/bBw+eUISCe3bSc/bB7ZNfZ32Otx4vNOslz0nPHS+IH4wLMb5nzlRKV8lO9R8H40YVcP9tF6V0HaDCITJgpUX5lBblc+ak0oTPd3Y6ew4d/WCLZdu+Ng4faaeru9J1CJXj3W6f+KT3GNf1VKL36HrB8ee8x5jE70XP8d3eu+e4E9fb47k+fp7u78VJz3mCn+fETKmUtiPZ0rCi0qL0/EpX4RBJkZwco6K0kIrSQmZVlUcdR2TIqJMnIiJJUeEQEZGkqHCIiEhSVDhERCQpKhwiIpIUFQ4REUmKCoeIiCRFhUNERJIyIq8AaGbNwHsDfPkEYPcQxhkqypUc5UqOciVnJOaa5u4VYQaOyMIxGGbWGPbyiemkXMlRruQoV3KyPZemqkREJCkqHCIikhQVjpPdHXWAXihXcpQrOcqVnKzOpR6HiIgkRVscIiKSlKwtHGZ2pZm9YWZvmdk3EjxfaGZLgueXmdn0DMn1e2bWbGarg+V/pSHTf5jZLjNb18vzZmbfDzKvMbO5qc4UMtdlZtbS7bP6VppyVZvZL81svZm9Zma3JxiT9s8sZK60f2ZmVmRmy83s1SDXXyYYk/bvY8hcaf8+dlt3rpm9YmaPJ3gutZ+Xu2fdAuQCbwOnAgXAq8BnmC0mAAAFV0lEQVTZPcZ8EbgruL0QWJIhuX4P+EGaP69LgbnAul6evwp4ivgVOC8ElmVIrsuAxyP4/2sKMDe4XQpsTPDfMe2fWchcaf/Mgs9gdHA7H1gGXNhjTBTfxzC50v597LbuPwMaEv33SvXnla1bHLXAW+6+yd2PAvcB1/QYcw3w0+D2g8DlluoLBofLlXbu/iKwt48h1wD/6XEvA+VmNiUDckXC3be7+6rg9gFgPTC1x7C0f2Yhc6Vd8BkcDO7mB0vP5mvav48hc0XCzKqATwM/7mVISj+vbC0cU4Et3e5v5eQv0Adj3L0daAHGZ0AugOuC6Y0Hzaw6xZnCCJs7ChcFUw1Pmdk56V55MEUwh/hfq91F+pn1kQsi+MyCaZfVwC5gqbv3+nml8fsYJhdE8338Z+BrQGcvz6f088rWwpGo8vb8SyLMmKEWZp2PAdPdfRbwHMf/qohSFJ9VGKuIn0bhPOBfgJ+nc+VmNhp4CPiKu+/v+XSCl6TlM+snVySfmbt3uPtsoAqoNbOZPYZE8nmFyJX276OZzQd2ufvKvoYleGzIPq9sLRxbge5/GVQBTb2NMbM8oIzUT4v0m8vd97j7keDuPcD5Kc4URpjPM+3cfX/XVIO7Pwnkm9mEdKzbzPKJ/3Je7O4/SzAkks+sv1xRfmbBOvcBLwBX9ngqiu9jv7ki+j5+GFhgZu8Sn87+uJn9d48xKf28srVwrADOMLNTzKyAePPo0R5jHgVuDm5fDzzvQacpylw95sEXEJ+njtqjwO8GewpdCLS4+/aoQ5nZ5K55XTOrJf7/+540rNeAfwfWu/t3exmW9s8sTK4oPjMzqzCz8uB2MTAP2NBjWNq/j2FyRfF9dPc73L3K3acT/x3xvLvf1GNYSj+vvKF6o+HE3dvN7I+BZ4jvyfQf7v6amf0foNHdHyX+BfsvM3uLeKVemCG5vmxmC4D2INfvpTqXmd1LfG+bCWa2Ffg28UYh7n4X8CTxvYTeAg4Dv5/qTCFzXQ/8kZm1A63AwjQUf4j/Rfg7wNpgfhzgz4FYt2xRfGZhckXxmU0BfmpmucQL1f3u/njU38eQudL+fexNOj8vHTkuIiJJydapKhERGSAVDhERSYoKh4iIJEWFQ0REkqLCISIiSVHhEOmDmf1P8O90M6sf4vf+80TrEsl02h1XJAQzuwz4qrvPT+I1ue7e0cfzB9199FDkE0knbXGI9MHMus6O+nfAR4JrLvxpcPK7O81sRXCCuy8E4y+z+DUvGoC1wWM/N7OVFr+mw63BY38HFAfvt7j7uoKjye80s3VmttbMbuj23i8EJ9PbYGaLu47yFkmnrDxyXGQAvkG3LY6gALS4+wVmVgj8xsyeDcbWAjPd/Z3g/i3uvjc4bcUKM3vI3b9hZn8cnECvp88Cs4HzgAnBa14MnpsDnEP8vFa/IX40+EtD/+OK9E5bHCID80ni55paTfzU5OOBM4LnlncrGhA/LcWrwMvETzx3Bn27BLg3ODPrTuBXwAXd3nuru3cCq4HpQ/LTiCRBWxwiA2PAn7j7Myc8GO+FHOpxfx5wkbsfNrMXgKIQ792bI91ud6DvsERAWxwi4RwgfrnVLs8QPxlgPoCZnWlmJQleVwa8HxSNGcQvE9vlWNfre3gRuCHoo1QQv0Tu8iH5KUSGgP5aEQlnDdAeTDn9P+B7xKeJVgUN6mbg2gSvexq4zczWAG8Qn67qcjewxsxWufuibo8/DFxE/JrzDnzN3XcEhUckctodV0REkqKpKhERSYoKh4iIJEWFQ0REkqLCISIiSVHhEBGRpKhwiIhIUlQ4REQkKSocIiKSlP8PWKaeHz5ON1sAAAAASUVORK5CYII=\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    },
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "-2.862398401331204\n"
     ]
    }
   ],
   "source": [
    "cost_val=[] # コスト関数の履歴\n",
    "init = np.random.rand(6)\n",
    "callback(init)\n",
    "res = scipy.optimize.minimize(cost, init,\n",
    "                              method='Powell',\n",
    "                              callback=callback)\n",
    "plt.plot(cost_val)\n",
    "plt.xlabel(\"iteration\")\n",
    "plt.ylabel(\"energy expectation value\")\n",
    "plt.show()\n",
    "print(cost(res.x))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "ノイズが小さければほとんど同じような振る舞いで最適化が行えることがわかる。（興味のある読者はぜひノイズを大きくして実験してみてほしい。）"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### 参考文献 \n",
    "[1] A. Peruzzo _et al_. , “A variational eigenvalue solver on a photonic quantum processor“ [Nat. Commun. 5:4213 doi: 10.1038/ncomms5213 (2014)](https://www.nature.com/articles/ncomms5213) "
   ]
  }
 ],
 "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.7.5"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}