select_r.ipynb

{
 "cells": [
  {
   "cell_type": "code",
   "execution_count": 1,
   "id": "814b334e",
   "metadata": {},
   "outputs": [],
   "source": [
    "from pyqpanda import *\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "144c67d8",
   "metadata": {},
   "source": [
    "## Init Quantum Environment"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "id": "a4420e15",
   "metadata": {},
   "outputs": [],
   "source": [
    "class InitQMachine:\n",
    "    def __init__(self, qubitsCount, cbitsCount, machineType = QMachineType.CPU):\n",
    "        self.machine = init_quantum_machine(machineType)\n",
    "        \n",
    "        self.qubits = self.machine.qAlloc_many(qubitsCount)\n",
    "        self.cbits = self.machine.cAlloc_many(cbitsCount)\n",
    "        \n",
    "        print(f'Init Quantum Machine with qubits:[{qubitsCount}] / cbits:[{cbitsCount}] Successfully')\n",
    "    \n",
    "    def __del__(self):\n",
    "        destroy_quantum_machine(self.machine)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "id": "4a6fee49",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Init Quantum Machine with qubits:[4] / cbits:[4] Successfully\n"
     ]
    }
   ],
   "source": [
    "ctx = InitQMachine(4, 4)\n",
    "\n",
    "machine = ctx.machine\n",
    "qubits = ctx.qubits\n",
    "cbits = ctx.cbits"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1440b9a3",
   "metadata": {},
   "source": [
    "# 0. Tool Functions"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "id": "31aa0408",
   "metadata": {},
   "outputs": [],
   "source": [
    "Dag = lambda matrix: matrix.conj().T"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "id": "6c6a5574",
   "metadata": {},
   "outputs": [],
   "source": [
    "normalize = lambda v: v / np.linalg.norm(v)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "id": "87a3febc",
   "metadata": {},
   "outputs": [],
   "source": [
    "def is_hermitian(matrix):\n",
    "    return np.allclose(matrix, Dag(matrix))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "id": "1ba9181c",
   "metadata": {},
   "outputs": [],
   "source": [
    "def ed(m, n):\n",
    "     return np.sqrt(np.sum((m - n) ** 2))"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "7b17c7db",
   "metadata": {},
   "source": [
    "# 1. Solve using numpy"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b3b505ab",
   "metadata": {},
   "source": [
    "## 1.0 solve original problem"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "id": "b89db3bf",
   "metadata": {},
   "outputs": [],
   "source": [
    "A = np.array([\n",
    "    [1, 1],\n",
    "    [2 ** 0.5 / 2, -(2 ** 0.5) / 2]\n",
    "])\n",
    "\n",
    "b = np.array([\n",
    "    [0.5], [-(2 ** 0.5) / 2]\n",
    "])"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "id": "17badd55",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.25],\n",
       "       [ 0.75]])"
      ]
     },
     "execution_count": 9,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# solve: Ax = b\n",
    "x = np.linalg.solve(A, b)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1c314ecd",
   "metadata": {},
   "source": [
    "## 1.1 test A†Ax = A†b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "id": "b63748ac",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[ 1.        ,  0.70710678],\n",
       "       [ 1.        , -0.70710678]])"
      ]
     },
     "execution_count": 10,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "Dag(A)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "id": "62b8296a",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_ = Dag(A) @ A # make A hermitian\n",
    "b_ = Dag(A) @ b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "id": "e65edf2c",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([[-0.25],\n",
       "       [ 0.75]])"
      ]
     },
     "execution_count": 12,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "x = np.linalg.solve(A_, b_)\n",
    "x"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "de5e76d7",
   "metadata": {},
   "source": [
    "## 1.2 eigenvalue of A"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "id": "9f31bed1",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "array([2., 1.])"
      ]
     },
     "execution_count": 13,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "eigenvalue, _ = np.linalg.eig(A_)\n",
    "eigenvalue"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8700682d",
   "metadata": {},
   "source": [
    "# 2. HHL"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f74258be",
   "metadata": {},
   "source": [
    "## 2.1 make A hermitian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "id": "8f95640a",
   "metadata": {},
   "outputs": [],
   "source": [
    "A_ = Dag(A) @ A # make A hermitian"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "id": "62c5d4ce",
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "False\n",
      "True\n"
     ]
    }
   ],
   "source": [
    "# test is_hermitian\n",
    "print(is_hermitian(A)) # false\n",
    "print(is_hermitian(A_)) # true"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "f244bbee",
   "metadata": {},
   "source": [
    "## 2.2 HHL algorithm subroutines"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "1f8b261b",
   "metadata": {},
   "source": [
    "### 2.2.1 encode b"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "id": "158724f9",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "[0j, (1+0j), (1+0j), 0j]"
      ]
     },
     "execution_count": 16,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "QOperator(X(qubits[0])).get_matrix()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "id": "2b3123b4",
   "metadata": {},
   "outputs": [],
   "source": [
    "def encode(b):\n",
    "    circuit = create_empty_circuit()\n",
    "#     circuit << amplitude_encode(qubits[3], b)\n",
    "    circuit << X(qubits[3])\n",
    "    \n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "id": "30dea6f4",
   "metadata": {},
   "outputs": [],
   "source": [
    "# https://arxiv.org/pdf/1110.2232.pdf"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "b8ea6301",
   "metadata": {},
   "source": [
    "### 2.2.2 phase estimation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "id": "e9d51043",
   "metadata": {},
   "outputs": [],
   "source": [
    "def phase_estimation(A):\n",
    "    circuit = create_empty_circuit()\n",
    "    \n",
    "    circuit << H(qubits[1]) << H(qubits[2]) << BARRIER(qubits[1:3])\n",
    "    \n",
    "#     circuit << QOracle(qubits[3], expMat(1j, A, np.pi / 2)).control(qubits[2]) # C-U^1\n",
    "#     circuit << QOracle(qubits[3], expMat(1j, A, np.pi)).control(qubits[1]) # C-U^2\n",
    "    circuit << CU(-np.pi / 4, -3 * np.pi / 2, -3 * np.pi / 2, 3 * np.pi / 2, qubits[2], qubits[3])\n",
    "    circuit << CU(-3 * np.pi/2, -3 * np.pi, -3 * np.pi, -2 * np.pi, qubits[1], qubits[3])\n",
    "    circuit << BARRIER(qubits[1:3])\n",
    "    \n",
    "    # inverse QFT\n",
    "    circuit << SWAP(qubits[1], qubits[2])\n",
    "    circuit << H(qubits[2])\n",
    "    circuit << S(qubits[2]).dagger().control(qubits[1])\n",
    "    circuit << H(qubits[1])\n",
    "    circuit << SWAP(qubits[1], qubits[2])\n",
    "\n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "a336901a",
   "metadata": {},
   "source": [
    "### 2.2.4 controlled rotations"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "id": "84ae536d",
   "metadata": {},
   "outputs": [],
   "source": [
    "def rotation(r):\n",
    "    circuit = create_empty_circuit()\n",
    "    \n",
    "    circuit << RY(qubits[0], 2*np.pi/(2**r)).control(qubits[1])\n",
    "    circuit << RY(qubits[0], np.pi/(2**r)).control(qubits[2])\n",
    "    \n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "8715c9d2",
   "metadata": {},
   "source": [
    "### 2.2.5 uncompute"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "id": "43b72335",
   "metadata": {},
   "outputs": [],
   "source": [
    "def uncompute(A):\n",
    "    circuit = create_empty_circuit()\n",
    "    \n",
    "    # QFT\n",
    "    circuit << SWAP(qubits[1], qubits[2])\n",
    "    circuit << H(qubits[1])\n",
    "    circuit << S(qubits[2]).control(qubits[1])\n",
    "    circuit << H(qubits[2])\n",
    "    circuit << SWAP(qubits[1], qubits[2])\n",
    "    circuit << BARRIER(qubits[1:3])\n",
    "\n",
    "#     circuit << QOracle(qubits[3], expMat(-1j, A, np.pi)).control(qubits[1])\n",
    "#     circuit << QOracle(qubits[3], expMat(-1j, A, np.pi / 2)).control(qubits[2])\n",
    "    circuit << CU(-3 * np.pi/2, -3 * np.pi, -3 * np.pi, -2 * np.pi, qubits[1], qubits[3])\n",
    "    circuit << CU(np.pi/4, -3*np.pi/2, -np.pi/2, -np.pi/2, qubits[2], qubits[3]) << BARRIER(qubits[1:3])\n",
    "    \n",
    "    circuit << H(qubits[1]) << H(qubits[2])\n",
    "    \n",
    "    return circuit"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "08711eb1",
   "metadata": {},
   "source": [
    "## 2.3 full HHL algorithm "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "id": "c6cff82b",
   "metadata": {},
   "outputs": [],
   "source": [
    "def HHL(A, b, r, flag=False):\n",
    "    prog = create_empty_qprog()\n",
    "    \n",
    "    # Step 0. check input\n",
    "    if not is_hermitian(A):\n",
    "        b = (Dag(A) @ b).round(4)\n",
    "        A = (Dag(A) @ A).round(4) # make A hermitian\n",
    "    \n",
    "    normed_b = (b / np.linalg.norm(b)).round(4)\n",
    "    \n",
    "    # Step 1. state preparation\n",
    "    prog << encode(normed_b)\n",
    "    \n",
    "    # Step 2. phase estimation\n",
    "    prog << phase_estimation(A)\n",
    "    \n",
    "    # Step 3. rotation\n",
    "    prog << rotation(r)\n",
    "    \n",
    "    # Step 4. uncompute\n",
    "    prog << uncompute(A)\n",
    "    \n",
    "    # Step 5. measure ancilla qubit\n",
    "    prog << Measure(qubits[0], cbits[0])\n",
    "    if flag:\n",
    "        results = [0, 0]\n",
    "        for i in range(10000):\n",
    "            result = prob_run_list(prog, qubits[0], -1)\n",
    "            results[0] += result[0]\n",
    "            results[1] += result[1]\n",
    "\n",
    "        results[0] /= 10000\n",
    "        results[1] /= 10000\n",
    "        \n",
    "        return results\n",
    "    \n",
    "    result = directly_run(prog)\n",
    "    \n",
    "    if not result['c0']:\n",
    "        return HHL(A, b, r)\n",
    "    \n",
    "    # Step 6. get results\n",
    "    prog << Measure(qubits[3], cbits[3])\n",
    "    \n",
    "    qstate = get_qstate()\n",
    "    normed_x = np.real(np.array([qstate[1], qstate[9]])) # 0001 1001\n",
    "        \n",
    "    # Step 7. recover x\n",
    "    N = len(normed_b)\n",
    "    ratio = 0.0\n",
    "    for i in range(N):\n",
    "        if not abs(normed_b[i]) < 1e-8:\n",
    "            ratio = normed_b[i][0] / np.sum([ normed_x[j] * A[i][j] for j in range(N) ])\n",
    "            break\n",
    "    \n",
    "    # normed_x = x / ||x|| => x = normed_x * ||x||\n",
    "    x_ = (normed_x * ratio)\n",
    "    distance = ed(x_, x.flatten())\n",
    "    \n",
    "    return distance, x_"
   ]
  },
  {
   "cell_type": "markdown",
   "id": "ae0abfae",
   "metadata": {},
   "source": [
    "## 2.4 select proper r for rotation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "id": "b4a014e5",
   "metadata": {},
   "outputs": [],
   "source": [
    "r = [1, 2, 3, 4, 5, 6, 7, 8]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "id": "26a25f1e",
   "metadata": {},
   "outputs": [],
   "source": [
    "p = [ HHL(A, b, i, True)[1] for i in r]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "id": "ed01e809",
   "metadata": {},
   "outputs": [],
   "source": [
    "distances = [ HHL(A, b, i)[0] for i in r ]"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "id": "141826ca",
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "id": "21c1b5cd",
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAABNcAAAF5CAYAAAC8xcwdAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABXPUlEQVR4nO3dfVxUdfr/8fcMCmjJpBI3Koaa5RIKioJYphWGrWtZu0WupWHr7rbaany3VdrS7EYszcVv+pN0syzXtCxNzcgiLV3pS4JsmWVq3pAKatqMYaHOzO8P1rFZQGFm4Azyej4e56Hzmeucc82Zdru65pzPx+R0Op0CAAAAAAAAUGdmoxMAAAAAAAAAGiuaawAAAAAAAICHaK4BAAAAAAAAHqK5BgAAAAAAAHiI5hoAAAAAAADgIZprAAAAAAAAgIdorgEAAAAAAAAeorkGAAAAAAAAeKiZ0Qn4C4fDoYMHD6pVq1YymUxGpwMAABoJp9OpEydOqF27djKb+d3SH1HnAQAAT9S2zqO59h8HDx5UVFSU0WkAAIBGqqSkRB06dDA6DVSDOg8AAHjjQnUezbX/aNWqlaTKCxYSEmJwNgAAoLGw2WyKiopy1RKQ5s6dqxkzZqi0tFRxcXF6/vnnlZiYWG3sF198ocmTJ6uwsFD79u3T3//+d02YMKHGY0+fPl2ZmZkaP368srOza5UPdR4AAPBEbes8mmv/cfYRgZCQEIouAABQZzxuWGnZsmXKyMhQTk6OkpKSlJ2drdTUVO3YsUNhYWFV4k+ePKnOnTvrzjvv1EMPPXTeY3/66ad64YUX1KNHjzrlRJ0HAAC8caE6j4lBAAAA4DOzZs3SmDFjlJ6erpiYGOXk5Khly5ZauHBhtfF9+vTRjBkzdPfddysoKKjG4/7www8aMWKEFixYoNatW9dX+gAAAHVGcw0AAAA+cerUKRUWFiolJcU1ZjablZKSovz8fK+OPXbsWA0ZMsTt2DWpqKiQzWZz2wAAAOoLzTUAAAD4xNGjR2W32xUeHu42Hh4ertLSUo+Pu3TpUhUVFSkrK6tW8VlZWbJYLK6NxQwAAEB9orkGAAAAv1VSUqLx48frn//8p4KDg2u1T2ZmpqxWq2srKSmp5ywBAEBTxoIGAAAA8InQ0FAFBASorKzMbbysrEwREREeHbOwsFCHDx9Wr169XGN2u10ff/yx5syZo4qKCgUEBLjtExQUdN752wAAAHyJO9cAAADgE4GBgUpISFBeXp5rzOFwKC8vT8nJyR4d86abbtLnn3+u4uJi19a7d2+NGDFCxcXFVRprAAAADY071xqA3eFUwZ5jOnziJ4W1ClZipzYKMJ9/GVcAAIDGKCMjQ6NGjVLv3r2VmJio7OxslZeXKz09XZI0cuRItW/f3jV/2qlTp7R9+3bX3w8cOKDi4mJdeumluvLKK9WqVSvFxsa6neOSSy5R27Ztq4wbwmGXdm6UrIckS6TUtb9kpuEHAEBTQnOtnuVuO6Spq7frkPUn11ikJVhThsZocGykgZkBAAD4Xlpamo4cOaLJkyertLRU8fHxys3NdS1ysH//fpnN5x6eOHjwoHr27Ol6PXPmTM2cOVMDBgzQhg0bGjr9uil6S3p9vHT823NjrTtId82Wet1hXF4AAKBBmZxOp9PoJPyBzWaTxWKR1WpVSEiIT46Zu+2QHlhcpP++wGfvWZt3Ty8abAAANHL1UUPAt+rlOyp6S3rhN1JNld4fltNgAwCgkattDcGca/XE7nBq6urtVcot6VwJNnX1dtkd9DYBAAAaFYe98o6181V6r0+ojAMAABc9mmv1pGDPMbdHQf+bU9Ih608q2HOs4ZICAACA93ZudH8UtAqndLykMg4AAFz0aK7Vk8Mnam6seRIHAAAAP2E95Ns4AADQqNFcqydhrYJ9GgcAAAA/YanlnLm1jQMAAI0azbV6ktipjSItwa7FC/6bSZWrhiZ2atOQaQEAAMBbXftXrgp6vkqvdVRlHAAAuOjRXKsnAWaTpgyNkVS17Dr7esrQGAWYayrKAAAA4JfMAdJds//zooZK767syjgAAHDRo7lWjwbHRmrePb0UYXF/9DPCEqx59/TS4FgeFQAAAGiUet0h/WG51Lq9+3jrDpXjve4wJi8AANDgmhmdwMVucGykBsVEqGDPMR0+8ZPCWlU+CsodawAAAI1crzuk+NsqVwW1HqqcY61rf+5YAwCgiaG51gACzCYld2lrdBoAAADwNXOAdPVAo7MAAAAG4rFQAAAAAAAAwEM01wAAAAAAAAAP0VwDAAAAAAAAPERzDQAAAAAAAPAQzTUAAAAAAADAQzTXAAAAAAAAAA/RXAMAAAAAAAA8RHMNAAAAAAAA8JBfN9fmzp2r6OhoBQcHKykpSQUFBTXGDhw4UCaTqco2ZMiQBswYAAAAAAAATYnfNteWLVumjIwMTZkyRUVFRYqLi1NqaqoOHz5cbfxbb72lQ4cOubZt27YpICBAd955ZwNnDgAAAAAAgKbCb5trs2bN0pgxY5Senq6YmBjl5OSoZcuWWrhwYbXxbdq0UUREhGt7//331bJlS5prAAAAAAAAqDd+2Vw7deqUCgsLlZKS4hozm81KSUlRfn5+rY7x4osv6u6779Yll1xS7fsVFRWy2WxuGwAAAAAAAFAXftlcO3r0qOx2u8LDw93Gw8PDVVpaesH9CwoKtG3bNv3ud7+rMSYrK0sWi8W1RUVFeZ03AAAAAAAAmha/bK5568UXX1T37t2VmJhYY0xmZqasVqtrKykpacAMAQAAAAAAcDFoZnQC1QkNDVVAQIDKysrcxsvKyhQREXHefcvLy7V06VI98cQT540LCgpSUFCQ17kCAAAAAACg6fLLO9cCAwOVkJCgvLw815jD4VBeXp6Sk5PPu+8bb7yhiooK3XPPPfWdJgAAAAAAAJo4v7xzTZIyMjI0atQo9e7dW4mJicrOzlZ5ebnS09MlSSNHjlT79u2VlZXltt+LL76oYcOGqW3btkakDQAAAAAAgCbEb5traWlpOnLkiCZPnqzS0lLFx8crNzfXtcjB/v37ZTa733i3Y8cObdq0SevWrTMiZQAAAAAAADQxJqfT6TQ6CX9gs9lksVhktVoVEhJidDoAAKCRoIbwf3xHAADAE7WtIfxyzjUAAAAAAACgMaC5BgAAAAAAAHiI5hoAAAB8au7cuYqOjlZwcLCSkpJUUFBQY+wXX3yhX//614qOjpbJZFJ2dnaVmKysLPXp00etWrVSWFiYhg0bph07dtTjJwAAAKg9mmsAAADwmWXLlikjI0NTpkxRUVGR4uLilJqaqsOHD1cbf/LkSXXu3FnTp09XREREtTEfffSRxo4dq08++UTvv/++Tp8+rZtvvlnl5eX1+VEAAABqhQUN/oOJbgEAgCeoIdwlJSWpT58+mjNnjiTJ4XAoKipKDz74oCZNmnTefaOjozVhwgRNmDDhvHFHjhxRWFiYPvroI11//fUXzInvCAAAeIIFDQAAANCgTp06pcLCQqWkpLjGzGazUlJSlJ+f77PzWK1WSVKbNm2qfb+iokI2m81tAwAAqC801wAAAOATR48eld1uV3h4uNt4eHi4SktLfXIOh8OhCRMm6Nprr1VsbGy1MVlZWbJYLK4tKirKJ+cGAACoDs01AAAANBpjx47Vtm3btHTp0hpjMjMzZbVaXVtJSUkDZggAAJqaZkYnAAAAgItDaGioAgICVFZW5jZeVlZW42IFdTFu3DitWbNGH3/8sTp06FBjXFBQkIKCgrw+HwAAQG1w5xoAAAB8IjAwUAkJCcrLy3ONORwO5eXlKTk52ePjOp1OjRs3TitWrNCHH36oTp06+SJdAAAAn+DONQAAAPhMRkaGRo0apd69eysxMVHZ2dkqLy9Xenq6JGnkyJFq3769srKyJFUugrB9+3bX3w8cOKDi4mJdeumluvLKKyVVPgq6ZMkSvf3222rVqpVr/jaLxaIWLVoY8CkBAADOobkGAAAAn0lLS9ORI0c0efJklZaWKj4+Xrm5ua5FDvbv3y+z+dzDEwcPHlTPnj1dr2fOnKmZM2dqwIAB2rBhgyRp3rx5kqSBAwe6neull17SfffdV6+fBwAA4EJMTqfTaXQS/sBms8lischqtSokJMTodAAAQCNBDeH/+I4AAIAnaltDMOcaAAAAAAAA4CGaawAAAAAAAICHaK4BAAAAAAAAHqK5BgAAAAAAAHiI5hoAAAAAAADgIZprAAAAAAAAgIdorgEAAAAAAAAeorkGAAAAAAAAeIjmGgAAAAAAAOAhmmsAAAAAAACAh2iuAQAAAAAAAB6iuQYAAAAAAAB4iOYaAAAAAAAA4CG/ba7NnTtX0dHRCg4OVlJSkgoKCs4b//3332vs2LGKjIxUUFCQrrrqKq1du7aBsgUAAAAAAEBT1MzoBKqzbNkyZWRkKCcnR0lJScrOzlZqaqp27NihsLCwKvGnTp3SoEGDFBYWpuXLl6t9+/bat2+fLrvssoZPHgAAAAAAAE2GXzbXZs2apTFjxig9PV2SlJOTo3feeUcLFy7UpEmTqsQvXLhQx44d0+bNm9W8eXNJUnR0dEOmDAAAAAAAgCbI7x4LPXXqlAoLC5WSkuIaM5vNSklJUX5+frX7rFq1SsnJyRo7dqzCw8MVGxuradOmyW6313ieiooK2Ww2tw0AAAAAAACoC79rrh09elR2u13h4eFu4+Hh4SotLa12n2+++UbLly+X3W7X2rVr9dhjj+m5557TU089VeN5srKyZLFYXFtUVJRPPwcAAAAAAAAufn7XXPOEw+FQWFiY5s+fr4SEBKWlpelvf/ubcnJyatwnMzNTVqvVtZWUlDRgxgAAAAAAALgY+N2ca6GhoQoICFBZWZnbeFlZmSIiIqrdJzIyUs2bN1dAQIBr7Be/+IVKS0t16tQpBQYGVtknKChIQUFBvk0eAAAAAAAATYrf3bkWGBiohIQE5eXlucYcDofy8vKUnJxc7T7XXnutdu3aJYfD4Rr7+uuvFRkZWW1jDQAAAAAAAPAFv2uuSVJGRoYWLFigRYsW6csvv9QDDzyg8vJy1+qhI0eOVGZmpiv+gQce0LFjxzR+/Hh9/fXXeueddzRt2jSNHTvWqI8AAAAAAACAJsDvHguVpLS0NB05ckSTJ09WaWmp4uPjlZub61rkYP/+/TKbz/UFo6Ki9N577+mhhx5Sjx491L59e40fP14TJ0406iMAAAAAAACgCTA5nU6n0Un4A5vNJovFIqvVqpCQEKPTAQAAjQQ1hP/jOwIAAJ6obQ3hl4+FAgAAAAAAAI0BzTUAAAAAAADAQzTXAAAAAAAAAA/RXAMAAAAAAAA8RHMNAAAAAAAA8BDNNQAAAAAAAMBDNNcAAAAAAAAAD9FcAwAAgE/NnTtX0dHRCg4OVlJSkgoKCmqM/eKLL/TrX/9a0dHRMplMys7O9vqYAAAADYnmGgAAAHxm2bJlysjI0JQpU1RUVKS4uDilpqbq8OHD1cafPHlSnTt31vTp0xUREeGTYwIAADQkmmsAAADwmVmzZmnMmDFKT09XTEyMcnJy1LJlSy1cuLDa+D59+mjGjBm6++67FRQU5JNjVlRUyGazuW0AAAD1heYaAAAAfOLUqVMqLCxUSkqKa8xsNislJUX5+fkNdsysrCxZLBbXFhUV5dG5AQAAaoPmGgAAAHzi6NGjstvtCg8PdxsPDw9XaWlpgx0zMzNTVqvVtZWUlHh0bgAAgNpoZnQCAAAAgC8FBQXV+IgpAACAr3HnGgAAAHwiNDRUAQEBKisrcxsvKyurcbECI44JAADgSzTXAAAA4BOBgYFKSEhQXl6ea8zhcCgvL0/Jycl+c0wAAABf4rFQAAAA+ExGRoZGjRql3r17KzExUdnZ2SovL1d6erokaeTIkWrfvr2ysrIkVS5YsH37dtffDxw4oOLiYl166aW68sora3VMAAAAI9FcAwAAgM+kpaXpyJEjmjx5skpLSxUfH6/c3FzXggT79++X2Xzu4YmDBw+qZ8+ertczZ87UzJkzNWDAAG3YsKFWxwQAADCSyel0Oo1Owh/YbDZZLBZZrVaFhIQYnQ4AAGgkqCH8H98RAADwRG1rCOZcAwAAAAAAADxEcw0AAAAAAADwEM01AAAAAAAAwEM01wAAAAAAAAAP0VwDAAAAAAAAPERzDQAAAAAAAPAQzTUAAAAAAADAQzTXAAAAAAAAAA951VybOnWqvv32W1/lAgAAAAAAADQqXjfXOnXqpKFDh2rVqlVyOBy+ykuSNHfuXEVHRys4OFhJSUkqKCioMfbll1+WyWRy24KDg32aDwAAAAAAAPBzXjXXnnrqKXXs2FHvvPOObr/9dkVFRemxxx7T3r17vU5s2bJlysjI0JQpU1RUVKS4uDilpqbq8OHDNe4TEhKiQ4cOubZ9+/Z5nQcAAAAAAABQE6+aa4888oh2796tdevW6c4779R3332np59+WldeeaUGDx6sN998U2fOnPHo2LNmzdKYMWOUnp6umJgY5eTkqGXLllq4cGGN+5hMJkVERLi28PDwGmMrKipks9ncNgAAAAAAAKAufLKgQUpKipYuXaoDBw5o5syZuvrqq7Vu3Trddddd6tChgyZNmqSdO3fW+ninTp1SYWGhUlJSziVqNislJUX5+fk17vfDDz/oiiuuUFRUlG677TZ98cUXNcZmZWXJYrG4tqioqFrnBwAAAAAAAEg+Xi20bdu2ysjI0BdffKFNmzZp+PDhOnz4sGbMmKFu3brppptu0ooVKy54nKNHj8put1e58yw8PFylpaXV7nP11Vdr4cKFevvtt7V48WI5HA7169evxgUXMjMzZbVaXVtJSUndPzAAAAAAAACatGb1cdDdu3dr9erVysvLc4116NBB69ev14YNG5SQkKA333zTp3eLJScnKzk52fW6X79++sUvfqEXXnhBTz75ZJX4oKAgBQUF+ez8AAAAAAAAaHp8dufa6dOntXTpUt1000266qqr9Mwzz+jMmTPKyMjQV199pX379ulf//qXbrnlFm3ZskXjxo2r8VihoaEKCAhQWVmZ23hZWZkiIiJqlU/z5s3Vs2dP7dq1y6vPBQAAAAAAANTE6+bal19+qYyMDLVr104jRozQ+vXrlZycrFdeeUXffvutZs6cqa5du0qqvLtszZo1SkxM1EcffVTjMQMDA5WQkOB255vD4VBeXp7b3WnnY7fb9fnnnysyMtK7DwgAAAAAAADUwKvHQq+77jrl5+fL6XQqJCREDzzwgP74xz8qNjb2vPtdc801+vTTT88bk5GRoVGjRql3795KTExUdna2ysvLlZ6eLkkaOXKk2rdvr6ysLEnSE088ob59++rKK6/U999/rxkzZmjfvn363e9+581HBAAAAAAAAGrkVXNt8+bN6tWrl/74xz/qt7/9rVq2bFmr/X73u9/p+uuvP29MWlqajhw5osmTJ6u0tFTx8fHKzc11LXKwf/9+mc3nbrw7fvy4xowZo9LSUrVu3VoJCQnavHmzYmJiPP+AAAAAF7F169Zp3rx5Kigo0NGjR3XPPffoxRdflCS99957eu+99/SXv/xF7dq1MzhTAAAA/2VyOp1OT3cuLCxUQkKCL/MxjM1mk8VikdVqVUhIiNHpAACARqKx1hDjx4/XnDlz5HQ6demll+qHH37Qfffdp4ULF0qSPvvsM8XHx+u5557TQw89ZHC23mms3xEAADBWbWsIr+Zce+edd7Rq1aoLxq1evVpPPPGEN6cCAACAj7zyyit6/vnnlZCQoKKiItlstioxPXr0UFRUlFavXm1AhgAAAI2HV821xx9/XCtXrrxg3KpVqzR16lRvTgUAAAAfmTdvni677DK98847io+PrzGuR48e+uabbxouMQAAgEbI69VCa8Nut7vNjwYAAADjbNu2Tf369dPll19+3jiLxaKysrIGygoAAKBxapCO1xdffKHWrVs3xKkAAABQCyaT6YIxBw8eVIsWLRogGwAAgMarzquFjh492u31pk2bqoyddebMGe3YsUNbtmzRsGHDPEoQAAAAvtW1a1cVFRXp9OnTat68ebUxJ06cUHFxsa655poGzg4AAKBxqXNz7eWXX3b93WQyadeuXdq1a9d59+nRo4dmzJhR5+QAAADge3feeaf+9re/adKkSXruueeqjcnMzJTVatXdd9/dwNkBAAA0LnVurq1fv16S5HQ6deONN2rw4MGaOHFitbGBgYFq166drrjiCu+yBAAAgM9MmDBBS5cuVXZ2tjZv3qzbbrtNkrR79279/e9/14oVK7Rp0yb16tVLY8aMMThbAAAA/2ZyOp1OT3dOT09X//79a3wstDGx2WyyWCyyWq0KCQkxOh0AANBINNYa4siRI7rvvvv07rvvymQy6b9LwkGDBmnx4sUXXPSgMWis3xEAADBWbWuIOt+59nMvvfSSN7sDAADAIJdffrneeecd/fvf/9a6deu0d+9eORwOdejQQYMGDVJiYqLRKQIAADQKXjXXAAAA0LjFxcUpLi7O6DQAAAAaLXNdgjt37qwuXbpoz549rte13bp06VIvHwAAAAD+Ze7cuYqOjlZwcLCSkpJUUFBw3vg33nhD3bp1U3BwsLp37661a9e6vf/DDz9o3Lhx6tChg1q0aKGYmBjl5OTU50cAAACotTo11/bu3as9e/bo9OnTrte13c425AAAAGCsf/zjH2rTpo1yc3NrjHn33XfVpk0bt5Xia2PZsmXKyMjQlClTVFRUpLi4OKWmpurw4cPVxm/evFnDhw/X/fffr61bt2rYsGEaNmyYtm3b5orJyMhQbm6uFi9erC+//FITJkzQuHHjtGrVqjrlBgAAUB+8WtDgYsJEtwAAwBONsYZISUnRF198oQMHDshsrv63Vrvdrvbt26tHjx5at25drY+dlJSkPn36aM6cOZIkh8OhqKgoPfjgg5o0aVKV+LS0NJWXl2vNmjWusb59+yo+Pt51d1psbKzS0tL02GOPuWISEhJ0yy236KmnnrpgTo3xOwIAAMarbQ1RpzvXAAAA0Pht375dPXr0qLGxJkkBAQGKi4vT9u3ba33cU6dOqbCwUCkpKa4xs9mslJQU5efnV7tPfn6+W7wkpaamusX369dPq1at0oEDB+R0OrV+/Xp9/fXXuvnmm6s9ZkVFhWw2m9sGAABQX2iuAQAANDHHjh1TaGjoBeNCQ0N19OjRWh/36NGjstvtCg8PdxsPDw9XaWlptfuUlpZeMP75559XTEyMOnTooMDAQA0ePFhz587V9ddfX+0xs7KyZLFYXFtUVFStPwMAAEBd0VwDAABoYkJDQ7Vz584Lxu3cuVOtW7dugIzO7/nnn9cnn3yiVatWqbCwUM8995zGjh2rDz74oNr4zMxMWa1W11ZSUtLAGQMAgKakWV2CAwICPD6RyWTSmTNnPN4fAAAAvnHdddfpjTfe0IYNGzRw4MBqYzZs2KAtW7bojjvuqPVxQ0NDFRAQoLKyMrfxsrIyRUREVLtPRETEeeN//PFHPfLII1qxYoWGDBkiSerRo4eKi4s1c+bMKo+USlJQUJCCgoJqnTcAAIA36nTnWlRUlDp27OjRxu34AAAA/iEjI0Mmk0nDhg3TzJkzZbVaXe/ZbDbNnDlTt99+u8xmsx566KFaHzcwMFAJCQnKy8tzjTkcDuXl5Sk5ObnafZKTk93iJen99993xZ8+fVqnT5+uMj9cQECAHA5HrXMDAACoL3W6c23v3r31lAYAAAAaSmJiop577jllZGRo4sSJmjhxotq0aSOpcj62s2bMmKFrr722TsfOyMjQqFGj1Lt3byUmJio7O1vl5eVKT0+XJI0cOVLt27dXVlaWJGn8+PEaMGCAnnvuOQ0ZMkRLly7Vli1bNH/+fElSSEiIBgwYoIcfflgtWrTQFVdcoY8++kivvPKKZs2a5YvLAQAA4JU6NdcAAABwcRg/frx69uyp6dOn66OPPtJ3330nSWrRooUGDhyoiRMn1rhgwPmkpaXpyJEjmjx5skpLSxUfH6/c3FzXogX79+93uwutX79+WrJkiR599FE98sgj6tq1q1auXKnY2FhXzNKlS5WZmakRI0bo2LFjuuKKK/T000/rj3/8o5dXAQAAwHsmp9PpNDoJf2Cz2WSxWGS1WhUSEmJ0OgAAoJG4GGoIu93uaq6FhoZWeQSzsbsYviMAANDwaltD1OnOtf3790uS2rdvr4CAANfr2urYsWOd4gEAAFD/AgICFBYWZnQaAAAAjVKdmmvR0dEym83avn27rrrqKkVHR8tkMtVqX1YLBQAAAAAAwMWmTs2166+/XiaTSS1btnR7DQAAgMZl+/btevbZZ/Xxxx/r0KFDOnXqVLVx/EAKAABwfnVqrm3YsOG8rwEAAOD/8vPzlZKSoh9//FGS1KZNG0VERBicFQAAQOPEaqEAAABNTGZmpn788UdNmDBBjz76qNq0aWN0SgAAAI2Wz5eCOn78uI4fPy5fLEI6d+5cRUdHKzg4WElJSSooKKjVfkuXLpXJZNKwYcO8zgEAAOBis2XLFsXHx2vWrFk01gAAALzkk+baqlWrdPPNN+vSSy9VaGioQkND1apVK9188816++23PTrmsmXLlJGRoSlTpqioqEhxcXFKTU3V4cOHz7vf3r179Ze//EX9+/f36LwAAAAXu8DAQHXr1s3oNAAAAC4KXjXXnE6nRo8erdtvv10ffPCBTp48KYvFIovFopMnT+qDDz7QHXfcofvuu6/Od7LNmjVLY8aMUXp6umJiYpSTk6OWLVtq4cKFNe5jt9s1YsQITZ06VZ07d/bmowEAAFy0rrvuOm3fvt3oNAAAAC4KXjXXZs+erZdfflmRkZGaN2+evv/+ex07dkzHjh2T1WpVTk6OIiMj9eqrr2r27Nm1Pu6pU6dUWFiolJSUc4mazUpJSVF+fn6N+z3xxBMKCwvT/ffff8FzVFRUyGazuW0AAABNwbRp07Rr1y7NnTvX6FQAAAAaPa8WNJg/f75atmypjRs3qlOnTm7vtWrVSr///e81aNAgde/eXfPnz9eECRNqddyjR4/KbrcrPDzcbTw8PFxfffVVtfts2rRJL774ooqLi2t1jqysLE2dOrVWsQAAABeToqIipaen689//rNef/11DRo0SB06dJDZXP3vriNHjmzgDAEAABoPr5pre/bs0c0331ylsfZznTp10k033aR169Z5c6rzOnHihO69914tWLBAoaGhtdonMzNTGRkZrtc2m01RUVH1lSIAAIDfuO+++2QymeR0OrVx40Zt3LhRJpOpSpzT6ZTJZKK5BgAAcB5eNdcuv/xyBQYGXjCuefPmtW56SVJoaKgCAgJUVlbmNl5WVqaIiIgq8bt379bevXs1dOhQ15jD4ZAkNWvWTDt27FCXLl3c9gkKClJQUFCtcwIAALhYTJ48udpmGgAAAOrOq+ba7bffrsWLF+v48eNq3bp1tTHHjh3Thx9+qBEjRtT6uIGBgUpISFBeXp6GDRsmqbJZlpeXp3HjxlWJ79atmz7//HO3sUcffVQnTpzQ7NmzuSMNAADgZx5//HGjUwAAALhoeLWgwVNPPaXOnTvrxhtv1Icffljl/fXr12vQoEHq0qWLpk2bVqdjZ2RkaMGCBVq0aJG+/PJLPfDAAyovL1d6erqkyrk/MjMzJUnBwcGKjY112y677DK1atVKsbGxtbq7DgAAAAAAAKirOt25duONN1YZCwwMVGFhoQYNGqQ2bdroiiuukCTt379f3333nSSpb9++GjZsmPLy8mp9rrS0NB05ckSTJ09WaWmp4uPjlZub61rkYP/+/TVOugsAAAAAAAA0BJPT6XTWNtibZpbJZJLdbvd4//pms9lksVhktVoVEhJidDoAAKCRaMw1xKZNm/T2229r586dOnHihKorC00mU51+IPVHjfk7AgAAxqltDVGnO9f27NnjdWIAAAAwltPp1P33369Fixa5GmpnVw896+xrFj4AAAA4vzo1184+8gkAAIDGKycnRy+//LJ69+6trKwszZs3TytWrNCOHTv0zTffaNmyZXr11VeVkZGhP/3pT0anCwAA4Ne8Wi0UAAAAjc/LL7+sSy65RO+++67atm2rxYsXS5K6du2qrl27KjU1Vb/85S+Vlpamfv368QMrAADAefi0ufb999/XOF+HJHXs2NGXpwMAAIAHvvzyS/Xr109t27aVJNejn3a7XQEBAZKk3/zmN0pISNDMmTN12223GZYrAACAv/O6uVZaWqpHH31Uq1atcq0OWh2TyaQzZ854ezoAAAB4yeFwuBprktSyZUtJ0vHjxxUaGuoa79q1q955550Gzw8AAKAx8Xz5T0mHDh1S7969tXDhQgUFBenyyy+X0+lU3759FRYW5rqDLTk5Wf379/dJwgAAAPBO+/btdfDgQdfrs499bt261S3u66+/VrNmzCICAABwPl4115566ikdPHhQTzzxhEpKSnTLLbfIZDLpX//6lw4dOqQNGzaoW7duMplMevfdd32VMwAAALzQq1cvbd++XXa7XZJ08803y+l06q9//au++uornThxQjNmzFBhYaF69uxpcLYAAAD+zavmWm5urjp16qRHH3202vevv/56rVu3Tlu3btWTTz7pzakAAADgI7feequOHj3qeuQzLi5Od999t/7973/rmmuu0WWXXaZJkyapWbNmevrppw3OFgAAwL951Vw7cOCA4uPjXa/PToBbUVHhGmvfvr1uuOEGvf76696cCgAAAD4yfPhw/fjjjxoyZIhrbNGiRZo2bZr69OmjK6+8Ur/85S+Vl5enxMREAzMFAADwf15NohESEuL2+rLLLpNU2XTr3Lmzazw4OFgHDhzw5lQAAADwoaCgILfXzZs316RJkzRp0iSDMgIAAGicvLpzrWPHjtq/f7/rdWxsrCRp7dq1rrGTJ0/qX//6lyIjI705FQAAAHzklVde0ebNmy8Y98knn+iVV15pgIwAAAAaL6+aazfeeKM+++wzHTlyRFLl/B2XXHKJHn74YU2aNEnPP/+8brjhBpWVlemWW27xScIAAADwzn333ad//OMfF4x78cUXlZ6e3gAZAQAANF5ePRY6YsQIlZSUaPv27RowYIDatGmjF154Qenp6Xr22WdlMpnkdDp1zTXXMBkuAABAI+NwOGQymYxOAwAAwK951VyLi4vTa6+95jY2fPhwXXvttVq7dq2OHz+uq666SrfeequaN2/uVaIAAABoWN98802VOXYBAADgzqvmWk06duyoP/7xj/VxaAAAAHjgiSeecHtdXFxcZeysM2fOaMeOHfr44481aNCgOp9r7ty5mjFjhkpLSxUXF6fnn3/+vKuOvvHGG3rssce0d+9ede3aVc8884x++ctfusV8+eWXmjhxoj766COdOXNGMTExevPNN9WxY8c65wcAAOBLPm+uHT9+XFLlyqE8RgAAAOAfHn/8cdeUHSaTScXFxSouLj7vPmFhYZo2bVqdzrNs2TJlZGQoJydHSUlJys7OVmpqqnbs2KGwsLAq8Zs3b9bw4cOVlZWlX/3qV1qyZImGDRumoqIi12JZu3fv1nXXXaf7779fU6dOVUhIiL744gsFBwfXKTcAAID6YHI6nU5vD7Jq1SrNmTNHmzdv1o8//ihJatGihfr166exY8fqtttu8zrR+maz2WSxWGS1Wnn8AQAA1FpjqSEWLVokSXI6nRo9erSrWVWdwMBAtWvXTn379lVQUFCdzpOUlKQ+ffpozpw5kirnbYuKitKDDz6oSZMmVYlPS0tTeXm51qxZ4xrr27ev4uPjlZOTI0m6++671bx5c7366qt1yuWsxvIdAQAA/1LbGsKrO9ecTqfuv/9+LVq0SGd7dJdddpkk6fvvv9cHH3ygvLw83XvvvXrppZe4kw0AAMAgo0aNcv190aJFuuWWW9zGfOHUqVMqLCxUZmama8xsNislJUX5+fnV7pOfn6+MjAy3sdTUVK1cuVJSZXPunXfe0V//+lelpqZq69at6tSpkzIzMzVs2LBqj1lRUaGKigrXa5vN5t0HAwAAOA+zNzvPnj1bL7/8siIjIzVv3jx9//33OnbsmI4dOyar1aqcnBxFRkbq1Vdf1ezZs32VMwAAALywfv16/fWvf/X5cY8ePSq73a7w8HC38fDwcJWWlla7T2lp6XnjDx8+rB9++EHTp0/X4MGDtW7dOt1+++2644479NFHH1V7zKysLFksFtcWFRXlg08HAABQPa+aa/Pnz1fLli21ceNG/eEPf3C7Ra5Vq1b6/e9/r40bN6pFixaaP3++18ni4mF3OJW/+zu9XXxA+bu/k93h9dPJAACglk6ePKn9+/ervLzcbfz48eOaNGmSfvWrX+lPf/qTdu/ebVCG5zgcDknSbbfdpoceekjx8fGuHM8+NvrfMjMzZbVaXVtJSUlDpgwAAJoYrx4L3bNnj26++WZ16tSpxphOnTrppptu0rp167w5FS4iudsOaerq7Tpk/ck1FmkJ1pShMRocG2lgZgAANA1PPvmknn32WRUUFCghIUFS5aOUffv21a5du1zTfSxfvlz//ve/FRlZu38/h4aGKiAgQGVlZW7jZWVlioiIqHafiIiI88aHhoaqWbNmiomJcYv5xS9+oU2bNlV7zKCgoDrPFQcAAOApr+5cu/zyyxUYGHjBuObNmys0NNSbU+EikbvtkB5YXOTWWJOkUutPemBxkXK3HTIoMwAAmo4PP/xQXbp0cTXWJGnx4sXauXOnbrjhBr333nv685//rKNHj+rvf/97rY8bGBiohIQE5eXlucYcDofy8vKUnJxc7T7Jyclu8ZL0/vvvu+IDAwPVp08f7dixwy3m66+/1hVXXFHr3AAAAOqLV82122+/XR9++KGOHz9eY8yxY8f04Ycf1jjhLJoOu8Opqau3q7oHQM+OTV29nUdEAQCoZ/v371fXrl3dxlatWiWTyaSXXnpJgwYNUnZ2tq666iq9++67dTp2RkaGFixYoEWLFunLL7/UAw88oPLycqWnp0uSRo4c6bbgwfjx45Wbm6vnnntOX331lR5//HFt2bJF48aNc8U8/PDDWrZsmRYsWKBdu3Zpzpw5Wr16tf70pz95cRUAAAB8w6vm2lNPPaXOnTvrxhtv1Icffljl/fXr12vQoEHq0qWLpk2b5s2pcBEo2HOsyh1rP+eUdMj6kwr2HGu4pAAAaIKOHz/uWuFdqlwBftOmTerRo4fb5P9xcXF1nq8sLS1NM2fO1OTJkxUfH6/i4mLl5ua6Fi3Yv3+/Dh06d6d6v379tGTJEs2fP19xcXFavny5Vq5cqdjYWFfM7bffrpycHD377LPq3r27/vGPf+jNN9/Udddd5+EVAAAA8J06zbl24403VhkLDAxUYWGhBg0apDZt2rhuz9+/f7++++47SVLfvn01bNiwKrf8o2k5fKLmxponcQAAwDMRERHas2eP63VhYaGOHz+ue++91y3OZDJ5dPxx48a53Xn2cxs2bKgyduedd+rOO+887zFHjx6t0aNHe5QPAABAfapTc626Yugsp9Op7777ztVQ+7n8/HyPizNcPMJaBfs0DgAAeCY+Pl5r1qzRypUrddNNN+nJJ5+UyWTSr371K7e4nTt3ql27dgZlCQAA0DjUqbn28184G8LcuXM1Y8YMlZaWKi4uTs8//7wSExOrjX3rrbc0bdo07dq1S6dPn1bXrl31P//zP1V+gYVxEju1UaQlWKXWn6qdd80kKcISrMRObRo6NQAAmpS//vWvWrNmjX79619LqvyRND4+3u0phbKyMv373//W8OHDjUoTAACgUahTc60hV2RatmyZMjIylJOTo6SkJGVnZys1NVU7duxQWFhYlfg2bdrob3/7m7p166bAwECtWbNG6enpCgsLU2pqaoPljZoFmE2aMjRGDywukklya7Cdva9xytAYBZi5yxEAgPrUr18/rVixQjNnztTRo0eVkJCgadOmyWw+Nx3va6+9platWmnw4MEGZgoAAOD/TE6n0y+XZkxKSlKfPn00Z84cSZXLuEdFRenBBx/UpEmTanWMXr16aciQIXryyScvGGuz2WSxWGS1WhUSEuJV7ji/3G2HNHX1drfFDSItwZoyNEaDYyMNzAwAgLqjhvB/fEcAAMATta0h6nTnWk3Kysq0cOFCbdy4UQcOHJAktW/fXtdff73S09Ndq0PV1qlTp1RYWOi2TLvZbFZKSory8/MvuL/T6dSHH36oHTt26Jlnnqk2pqKiQhUVFa7XNputTjnCc4NjIzUoJkIFe47p8ImfFNaq8lFQ7lgDAAAAAACNjdfNtTfffFOjR4/WDz/8oJ/fBPf555/rvffe0/Tp0/Xiiy+65vSojaNHj8put1dpyoWHh+urr76qcT+r1ar27duroqJCAQEB+n//7/9p0KBB1cZmZWVp6tSptc4JvhVgNim5S1uj0wAAAAAAAPCKV821LVu2aPjw4XI4HLr99tt17733Kjo6WiaTSXv37tWrr76qFStW6Le//a3+9a9/qXfv3r7Ku1qtWrVScXGxfvjhB+Xl5SkjI0OdO3fWwIEDq8RmZmYqIyPD9dpmsykqKqpe8wMAADBC586dZTKZ9MEHH6hTp07q3Llzrfc1mUzavXt3PWYHAADQuHnVXMvKypLdbtfy5ct1++23u73Xo0cP3XrrrVqxYoV+/etfa/r06Vq+fHmtjhsaGqqAgACVlZW5jZeVlSkiIqLG/cxms6688kpJlUvMf/nll8rKyqq2uRYUFKSgoKBa5QMAANCY7d27V5J0+vRpt9cAAADwnvnCITXbtGmT+vXrV6Wx9nO33367rr32Wm3cuLHWxw0MDFRCQoLy8vJcYw6HQ3l5eUpOTq71cRwOh9u8agAAAE2Rw+GQw+HQVVdd5fa6thsAAABq5tWda1arVR07drxgXMeOHfXpp5/W6dgZGRkaNWqUevfurcTERGVnZ6u8vFzp6emSpJEjR6p9+/bKysqSVHkXXe/evdWlSxdVVFRo7dq1evXVVzVv3ry6fzAAAAAAAACgFrxqrkVERGjr1q0XjCsuLj7v45zVSUtL05EjRzR58mSVlpYqPj5eubm5rkUO9u/fL7P53I135eXl+tOf/qRvv/1WLVq0ULdu3bR48WKlpaXV7UMBAAAAAAAAtWRy/nyJzzr6wx/+oH/84x+aOHGinnzySQUEBLi973Q69dhjjykrK0tjxoxRTk6O1wnXF5vNJovFIqvVqpCQEKPTAQAAjURjrCE2b96s9evX68svv9Tx48dlMpnUpk0bxcTE6IYbblBSUpLRKfpUY/yOAACA8WpbQ3jVXPv222/Vs2dPHTt2TB07dtRdd92l6OhoSdK+ffv0xhtvaO/evWrbtq2KiorUoUMHT09V7yi6AACAJxpTDfHZZ59p9OjRricP/rsMNJlMkqTExES9+OKLiomJafAc60Nj+o4AAID/aJDmmiR9/vnnGjFihLZt21Z5wP8UZWcP2717d/3zn/9UbGysN6epdxRdAADAE42lhvj000914403qry8XJdccoluueUWxcfHKzQ0VE6nU0ePHtXWrVv13nvvqby8XK1atdKGDRvUs2dPo1P3WmP5jgAAgH+pbQ3h1ZxrUmXz7LPPPtOGDRu0ceNGHTx4UJLUrl079e/fXwMHDvT2FAAAAPCC3W7XiBEjVF5ervvvv1/PPfdcjQWizWZTRkaGFi5cqN/+9rfavn2768dTAAAAVOVVc+2OO+5QZGSk5s6dq4EDB9JIAwAA8ENvv/22du3apbS0NC1YsOC8sSEhIfrHP/6hEydOaPny5Vq9erVuvfXWBsoUAACg8TFfOKRma9eu1XfffeerXAAAAFAPVq9eLbPZrGnTptV6n6ysLEnSypUr6ykrAACAi4NXzbVOnTqpvLzcV7kAAACgHhQWFurqq69Wp06dar1P586d1a1bNxUWFtZjZgAAAI2fV8214cOH66OPPlJpaamv8gEAAICPHTp0SFdddVWd97vqqqtc8+kCAACgel411zIzM9W/f38NGDBAK1as0OnTp32VFwAAAHzEarXKYrHUeb+QkBDZbLZ6yAgAAODi4dWCBldffbUcDodKSkr0m9/8RiaTSWFhYQoODq4SazKZtHv3bm9OBwAAAA+cOXNGZnPdf1M1m806c+ZMPWQEAABw8fCqubZ37163106nk0dEAQAAAAAA0GR41VxzOBy+ygMAAAD1aNGiRVq0aJHRaQAAAFx0vGquAQAAoHFwOp0e7WcymXycCQAAwMXFo+ba2rVrtXLlSpWUlCgoKEg9evRQenp6nZZ3BwAAQMPgaQMAAID6U+fm2ogRI7R06VJJ534BXb16tWbOnKmlS5fq1ltv9W2GAAAAAAAAgJ+qU3PtxRdf1GuvvaZmzZrp3nvvVc+ePXXixAmtWbNG+fn5GjlypPbt2+fRUu8AAAAAAABAY1On5tqiRYtkNpv17rvv6qabbnKNZ2ZmKj09Xa+88oreeustpaen+zxRAAAAAAAAwN+Y6xL8+eefq2/fvm6NtbMeeeQROZ1Off755z5LDgAAAAAAAPBndWqu2Ww2denSpdr3zo7bbDbvswIAAAAAAAAagTo115xOpwICAqo/kLnyUKxGBQAA0LTNnTtX0dHRCg4OVlJSkgoKCs4b/8Ybb6hbt24KDg5W9+7dtXbt2hpj//jHP8pkMik7O9vHWQMAAHimTs01AAAA4HyWLVumjIwMTZkyRUVFRYqLi1NqaqoOHz5cbfzmzZs1fPhw3X///dq6dauGDRumYcOGadu2bVViV6xYoU8++UTt2rWr748BAABQayan0+msbbDZbJbJZPLsRCaTzpw549G+DcFms8lischqtSokJMTodAAAQCNBDeEuKSlJffr00Zw5cyRVPtUQFRWlBx98UJMmTaoSn5aWpvLycq1Zs8Y11rdvX8XHxysnJ8c1duDAASUlJem9997TkCFDNGHCBE2YMKFWOfEdAQAAT9S2hqjznWtOp9OjjcdFAQAALm6nTp1SYWGhUlJSXGNms1kpKSnKz8+vdp/8/Hy3eElKTU11i3c4HLr33nv18MMP65prrrlgHhUVFbLZbG4bAABAfalTc83hcHi1AQAA4OJ19OhR2e12hYeHu42Hh4ertLS02n1KS0svGP/MM8+oWbNm+vOf/1yrPLKysmSxWFxbVFRUHT8JAABA7THnGgAAAPxWYWGhZs+erZdffrnW05NkZmbKarW6tpKSknrOEgAANGU01wAAAOAToaGhCggIUFlZmdt4WVmZIiIiqt0nIiLivPEbN27U4cOH1bFjRzVr1kzNmjXTvn379D//8z+Kjo6u9phBQUEKCQlx2wAAAOoLzTUAAAD4RGBgoBISEpSXl+caczgcysvLU3JycrX7JCcnu8VL0vvvv++Kv/fee/XZZ5+puLjYtbVr104PP/yw3nvvvfr7MAAAALXk1821uXPnKjo6WsHBwUpKSlJBQUGNsQsWLFD//v3VunVrtW7dWikpKeeNBwAAgO9lZGRowYIFWrRokb788ks98MADKi8vV3p6uiRp5MiRyszMdMWPHz9eubm5eu655/TVV1/p8ccf15YtWzRu3DhJUtu2bRUbG+u2NW/eXBEREbr66qsN+YwAAAA/57fNtWXLlikjI0NTpkxRUVGR4uLilJqaqsOHD1cbv2HDBg0fPlzr169Xfn6+oqKidPPNN+vAgQMNnDkAAEDTlZaWppkzZ2ry5MmKj49XcXGxcnNzXYsW7N+/X4cOHXLF9+vXT0uWLNH8+fMVFxen5cuXa+XKlYqNjTXqIwAAANSJyel0Oo1OojpJSUnq06eP5syZI6nykYKoqCg9+OCDmjRp0gX3t9vtat26tebMmaORI0deMN5ms8lischqtTIvBwAAqDVqCP/HdwQAADxR2xrCL+9cO3XqlAoLC5WSkuIaM5vNSklJUX5+fq2OcfLkSZ0+fVpt2rSp9v2KigrZbDa3DQAAAAAAAKgLv2yuHT16VHa73fX4wFnh4eEqLS2t1TEmTpyodu3auTXofi4rK0sWi8W1RUVFeZ03AAAAAAAAmha/bK55a/r06Vq6dKlWrFih4ODgamMyMzNltVpdW0lJSQNnCQAAAAAAgMaumdEJVCc0NFQBAQEqKytzGy8rK1NERMR59505c6amT5+uDz74QD169KgxLigoSEFBQT7JFwAAAAAAAE2TX965FhgYqISEBOXl5bnGHA6H8vLylJycXON+zz77rJ588knl5uaqd+/eDZEqAAAAAAAAmjC/vHNNkjIyMjRq1Cj17t1biYmJys7OVnl5udLT0yVJI0eOVPv27ZWVlSVJeuaZZzR58mQtWbJE0dHRrrnZLr30Ul166aWGfQ4AAAAAAABcvPy2uZaWlqYjR45o8uTJKi0tVXx8vHJzc12LHOzfv19m87kb7+bNm6dTp07pN7/5jdtxpkyZoscff7whUwcAAACaLodd2rlRsh6SLJFS1/6SOcDorAAAqDcmp9PpNDoJf2Cz2WSxWGS1WhUSEmJ0OgAAoJGghvB/fEcNqOgt6fXx0vFvz4217iDdNVvqdYdxeQEA4IHa1hB+OecaAAAAgEam6C3phd+4N9Yk6fiByvGit4zJCwCAekZzDQAAAIB3HPbKO9ZU3UMx/xl7fUJlHAAAFxmaawAAAAC8s3Nj1TvW3Dil4yWVcQAAXGRorgEAAADwjvWQb+MAAGhEaK4BAAAA8I4l0rdxAAA0IjTXAAAAAHina//KVUFlqiHAJLWOqowDAOAiQ3MNAAAAgHfMAdJds//z4r8bbP95fVd2ZRwAABcZmmsAAAAAvNfrDukPy6XW7d3HW3eoHO91hzF5AQBQz5oZnQAAAACAi0SvO6T42ypXBbUeqpxjrWt/7lgDAFzUaK4BAAAA8B1zgHT1QKOzAACgwfBYKAAAAAAAAOAhmmsAAAAAAACAh2iuAQAAAAAAAB6iuQYAAAAAAAB4iOYaAAAAAAAA4CGaawAAAAAAAICHaK4BAAAAAAAAHqK5BgAAAAAAAHiI5hoAAAAAAADgIZprAAAAAAAAgIdorgEAAAAAAAAeorkGAAAAAAAAeIjmGgAAAAAAAOAhmmsAAADwqblz5yo6OlrBwcFKSkpSQUHBeePfeOMNdevWTcHBwerevbvWrl3reu/06dOaOHGiunfvrksuuUTt2rXTyJEjdfDgwfr+GAAAALVCcw0AAAA+s2zZMmVkZGjKlCkqKipSXFycUlNTdfjw4WrjN2/erOHDh+v+++/X1q1bNWzYMA0bNkzbtm2TJJ08eVJFRUV67LHHVFRUpLfeeks7duzQrbfe2pAfCwAAoEYmp9PpNDoJf2Cz2WSxWGS1WhUSEmJ0OgAAoJGghnCXlJSkPn36aM6cOZIkh8OhqKgoPfjgg5o0aVKV+LS0NJWXl2vNmjWusb59+yo+Pl45OTnVnuPTTz9VYmKi9u3bp44dO1Z5v6KiQhUVFa7XNptNUVFRfEcAAKBOalvncecaAAAAfOLUqVMqLCxUSkqKa8xsNislJUX5+fnV7pOfn+8WL0mpqak1xkuS1WqVyWTSZZddVu37WVlZslgsri0qKqruHwYAAKCW/La5Vpe5Or744gv9+te/VnR0tEwmk7KzsxsuUQAAAEiSjh49KrvdrvDwcLfx8PBwlZaWVrtPaWlpneJ/+uknTZw4UcOHD6/xF+TMzExZrVbXVlJS4sGnAQAAqB2/bK7Vda6OkydPqnPnzpo+fboiIiIaOFsAAAA0hNOnT+uuu+6S0+nUvHnzaowLCgpSSEiI2wYAAFBf/LK5NmvWLI0ZM0bp6emKiYlRTk6OWrZsqYULF1Yb36dPH82YMUN33323goKCanWOiooK2Ww2tw0AAACeCw0NVUBAgMrKytzGy8rKavwBNCIiolbxZxtr+/bt0/vvv0/DDAAA+A2/a655MleHJ5iLA/7G7nAqf/d3erv4gPJ3fye7g7VGAACNS2BgoBISEpSXl+caczgcysvLU3JycrX7JCcnu8VL0vvvv+8Wf7axtnPnTn3wwQdq27Zt/XwAAAAADzQzOoH/dr65Or766iufnSczM1MZGRmu12dXkQKMkLvtkKau3q5D1p9cY5GWYE0ZGqPBsZEGZgYAQN1kZGRo1KhR6t27txITE5Wdna3y8nKlp6dLkkaOHKn27dsrKytLkjR+/HgNGDBAzz33nIYMGaKlS5dqy5Ytmj9/vqTKxtpvfvMbFRUVac2aNbLb7a752Nq0aaPAwEBjPigAAMB/+F1zraEEBQXV+hFSoD7lbjukBxYX6b/vUyu1/qQHFhdp3j29aLABABqNtLQ0HTlyRJMnT1Zpaani4+OVm5vr+uF0//79MpvPPTzRr18/LVmyRI8++qgeeeQRde3aVStXrlRsbKwk6cCBA1q1apUkKT4+3u1c69ev18CBAxvkcwEAANTE75prnszVATRWdodTU1dvr9JYkySnJJOkqau3a1BMhALMpgbODgAAz4wbN07jxo2r9r0NGzZUGbvzzjt15513VhsfHR0tp5OpEgAAgP/yuznXPJmrA2isCvYcc3sU9L85JR2y/qSCPccaLikAAAAAAFBrfnfnmlT3uTpOnTql7du3u/5+4MABFRcX69JLL9WVV15p2OcALuTwiZoba57EAQAAAACAhuWXzbW6ztVx8OBB9ezZ0/V65syZmjlzpgYMGFDtoweAvwhrFezTOAAAAAAA0LD8srkm1W2uDubiQGOV2KmNIi3BKrX+VO28ayZJEZZgJXZq09CpAQAAAACAWvC7OdeApiTAbNKUoTGSKhtpP3f29ZShMSxmAAAAAACAn6K5BhhscGyk5t3TSxEW90c/IyzBmndPLw2OjTQoMwAAAAAAcCF++1go0JQMjo3UoJgIFew5psMnflJYq8pHQbljDQAAAAAA/0ZzDfATAWaTkru0NToNAAAAAABQBzwWCgAAAAAAAHiI5hoAAAAAAADgIZprAAAAAAAAgIdorgEAAAAAAAAeorkGAAAAAAAAeIjVQgEAAACgPjns0s6NkvWQZImUuvaXzAFGZwUA8BGaawAAAABQX4rekl4fLx3/9txY6w7SXbOlXncYlxcAwGd4LBQAAAAA6kPRW9ILv3FvrEnS8QOV40VvGZMXAMCnaK4BAAAAgK857JV3rMlZzZv/GXt9QmUcAKBRo7kGAAAAAL62c2PVO9bcOKXjJZVxAIBGjeYaAAAAAPia9ZBv4wAAfovmGgAAAAD4miXSt3EAAL/FaqEA/JLd4VTBnmM6fOInhbUKVmKnNgowm4xOCwAAoHa69q9cFfT4AVU/75qp8v2u/Rs6MwCAj9FcA+B3crcd0tTV23XI+pNrLNISrClDYzQ4ll93AQBAI2AOkO6aXbkqqExyb7D95wfDu7Ir4wAAjRqPhQLwK7nbDumBxUVujTVJKrX+pAcWFyl3G/OSAACARqLXHdIflkut27uPt+5QOd7rDmPyAgD4FHeuAfAbdodTU1dvr3HBepOkqau3a1BMBI+IAgCAxqHXHVL8bZWrgloPVc6x1rU/d6wBwEWE5hoAv1Gw51iVO9Z+zinpkPUnFew5puQubRsuMQAAAG+YA6SrBxqdBQCgntBcA+A3Dp+oubHmSRwAAAD8kMPOnXwALio01wD4jbBWwT6NAwAAgJ8pekt6fbx0/NtzY607VC7+wBx0ABopFjQA4DcSO7VRpCVYNc2mZlLlqqGJndo0ZFp+w+5wKn/3d3q7+IDyd38nu6O62ekAAAD8VNFblaun/ryxJknHD1SOF71lTF4A4CXuXAPgNwLMJk0ZGqMHFhfVtGC9pgyNaZKLGeRuO6Spq7e7zUkXaQnWlKExGhwbaWBmxrE7nCrYc0yHT/yksFaVTdem+M8GAACNgsNeecfa+Zauen1C5eIPTfERUR6VrYprgkaE5hoAvzI4NlLz7ulVpZEU0YQbSbnbDumBxUVVStFS6096YHGR5t3Tq8ldF5qN1aPhWBXXBAD8xM6NVe9Yc+OUjpdUxjW1xR94VLYqrklVNBur8qNr4tePhc6dO1fR0dEKDg5WUlKSCgoKzhv/xhtvqFu3bgoODlb37t21du3aBsoUgC8Njo3Upok36rUxfTX77ni9NqavNk28sUk2TewOp6au3l7jb7ySNHX19ib1iOjZZuN/ryx7ttmYu+2QQZkZK3fbIV33zIcavuATjV9arOELPtF1z3zYZK+HxDUxkq9rOKfTqcmTJysyMlItWrRQSkqKdu7cWZ8fAYCvWWv5/721jbtY8KhsVVyTqorekh6JlmbdIL3428o/H4lumtfiLD+7Jn7bXFu2bJkyMjI0ZcoUFRUVKS4uTqmpqTp8+HC18Zs3b9bw4cN1//33a+vWrRo2bJiGDRumbdu2NXDmAHwhwGxScpe2ui2+vZK7tG2yd5oU7DlWpYn0c05Jh6w/qWDPsYZLykA0G6tHw7Eqrolx6qOGe/bZZ/W///u/ysnJ0f/93//pkksuUWpqqn76idWjgUbDUssfSWsbdzG44KOyqnxU1mFvwKQMxjWpimZjVX54Tfy2uTZr1iyNGTNG6enpiomJUU5Ojlq2bKmFCxdWGz979mwNHjxYDz/8sH7xi1/oySefVK9evTRnzpwGzhwAfOfwidr9h2Nt4xo7mo1V0XCsimtiLF/XcE6nU9nZ2Xr00Ud12223qUePHnrllVd08OBBrVy5sgE/GQCvdO1f+Vjf+Zauah1VGddU1OVR2aaCa+KOZmNVfnpN/LK5durUKRUWFiolJcU1ZjablZKSovz8/Gr3yc/Pd4uXpNTU1BrjKyoqZLPZ3DYA8DdhrYJ9GtfY0WysioZjVVwT49RHDbdnzx6Vlpa6xVgsFiUlJVHnAY2JOaByvixJVRts/3l9V3bTmkOKR2Wr4pq4o9lYlZ9eE79srh09elR2u13h4eFu4+Hh4SotLa12n9LS0jrFZ2VlyWKxuLaoqCjfJA8APpTYqY0iLcHn+41XkZbKSdqbApqNVdFwrIprYpz6qOHO/kmdB1wEet0h/WG51Lq9+3jrDpXjTW2ieh6VrYpr4o5mY1V+ek38srnWEDIzM2W1Wl1bSUmJ0SkBQBUBZpOmDI2RVONvvJoyNKbJzElHs7EqGo5VcU1AnQf4sV53SNP2ShnrpfuXVP45bU/Ta6xJPCpbHa6JO5qNVfnpNfHL5lpoaKgCAgJUVlbmNl5WVqaIiIhq94mIiKhTfFBQkEJCQtw2APBHg2MjNe+eXoqwuDcCIizBmndPrya1iirNxqpoOFbFNTFOfdRwZ/+kzgMuIuYA6eqBUuLwyj+b0qOgP8ejslVxTdzRbKzKT6+JXzbXAgMDlZCQoLy8PNeYw+FQXl6ekpOTq90nOTnZLV6S3n///RrjAaAxGRwbqU0Tb9RrY/pq9t3xem1MX22aeGOTaqydRbPRHQ3HqrgmxqmPGq5Tp06KiIhwi7HZbPq///s/6jwAjR+PylbFNTmHZmNVfnpNTE6n0y+Xylq2bJlGjRqlF154QYmJicrOztbrr7+ur776SuHh4Ro5cqTat2+vrKwsSZXLuA8YMEDTp0/XkCFDtHTpUk2bNk1FRUWKjY294PlsNpssFousViu/bgJAI2B3OFWw55gOn/hJYa0q70Jqys2S3G2HNHX1dreJ/CMtwZoyNKbJNRzPaqhrQg3hrj5quGeeeUbTp0/XokWL1KlTJz322GP67LPPtH37dgUHX/jxXr4jAH7PYa+cgN16qPJxtq79m1bDpDpck3OK3qpcIfPnE/m3jqpsIjWlZuPPNdA1qW0N0cxnZ/SxtLQ0HTlyRJMnT1Zpaani4+OVm5vrmsx2//79MpvP3XjXr18/LVmyRI8++qgeeeQRde3aVStXrqxVYw0A0PgEmE1K7tLW6DT8xuDYSA2KiaDh+DNcE2PURw3317/+VeXl5fr973+v77//Xtddd51yc3Nr1VgDgEbh7KOyOIdrck6vO6T422g2/pyfXRO/vXOtofGLJgAA8AQ1hP/jOwIAAJ6obQ3hl3OuAQAAAAAAAI0BzTUAAAAAAADAQzTXAAAAAAAAAA/RXAMAAAAAAAA8RHMNAAAAAAAA8BDNNQAAAAAAAMBDzYxOwF84nU5JlcusAgAA1NbZ2uFsLQH/Q50HAAA8Uds6j+baf5w4cUKSFBUVZXAmAACgMTpx4oQsFovRaaAa1HkAAMAbF6rzTE5+ZpUkORwOHTx4UK1atZLJZPL58W02m6KiolRSUqKQkBCfH78x4ppUxTVxx/WoimvijutRFdekqvq+Jk6nUydOnFC7du1kNjPjhj+izmt4XJOquCbuuB5VcU2q4pq443pU5S91Hneu/YfZbFaHDh3q/TwhISH8j+C/cE2q4pq443pUxTVxx/WoimtSVX1eE+5Y82/UecbhmlTFNXHH9aiKa1IV18Qd16Mqo+s8fl4FAAAAAAAAPERzDQAAAAAAAPAQzbUGEhQUpClTpigoKMjoVPwG16Qqrok7rkdVXBN3XI+quCZVcU1Q3/hnrCquSVVcE3dcj6q4JlVxTdxxParyl2vCggYAAAAAAACAh7hzDQAAAAAAAPAQzTUAAAAAAADAQzTXAAAAAAAAAA/RXAMAAAAAAAA8RHOtnn388ccaOnSo2rVrJ5PJpJUrVxqdkuGysrLUp08ftWrVSmFhYRo2bJh27NhhdFqGmTdvnnr06KGQkBCFhIQoOTlZ7777rtFp+Y3p06fLZDJpwoQJRqdimMcff1wmk8lt69atm9FpGe7AgQO655571LZtW7Vo0ULdu3fXli1bjE7LMNHR0VX+OTGZTBo7dqzRqRnCbrfrscceU6dOndSiRQt16dJFTz75pFjHCb5GreeOOs8ddd6FUetR61WHOs8ddZ47f6zzmhl25iaivLxccXFxGj16tO644w6j0/ELH330kcaOHas+ffrozJkzeuSRR3TzzTdr+/btuuSSS4xOr8F16NBB06dPV9euXeV0OrVo0SLddttt2rp1q6655hqj0zPUp59+qhdeeEE9evQwOhXDXXPNNfrggw9cr5s1a9r/9338+HFde+21uuGGG/Tuu+/q8ssv186dO9W6dWujUzPMp59+Krvd7nq9bds2DRo0SHfeeaeBWRnnmWee0bx587Ro0SJdc8012rJli9LT02WxWPTnP//Z6PRwEaHWc0ed54467/yo9c6h1juHOq8q6jx3/ljnNd3/xTaQW265RbfccovRafiV3Nxct9cvv/yywsLCVFhYqOuvv96grIwzdOhQt9dPP/205s2bp08++aRJF10//PCDRowYoQULFuipp54yOh3DNWvWTBEREUan4TeeeeYZRUVF6aWXXnKNderUycCMjHf55Ze7vZ4+fbq6dOmiAQMGGJSRsTZv3qzbbrtNQ4YMkVT5i+9rr72mgoICgzPDxYZazx11njvqvJpR67mj1juHOq8q6jx3/ljn8VgoDGe1WiVJbdq0MTgT49ntdi1dulTl5eVKTk42Oh1DjR07VkOGDFFKSorRqfiFnTt3ql27durcubNGjBih/fv3G52SoVatWqXevXvrzjvvVFhYmHr27KkFCxYYnZbfOHXqlBYvXqzRo0fLZDIZnY4h+vXrp7y8PH399deSpH//+9/atGkTTRCggVHnnUOd545azx213jnUeedHneefdR53rsFQDodDEyZM0LXXXqvY2Fij0zHM559/ruTkZP3000+69NJLtWLFCsXExBidlmGWLl2qoqIiffrpp0an4heSkpL08ssv6+qrr9ahQ4c0depU9e/fX9u2bVOrVq2MTs8Q33zzjebNm6eMjAw98sgj+vTTT/XnP/9ZgYGBGjVqlNHpGW7lypX6/vvvdd999xmdimEmTZokm82mbt26KSAgQHa7XU8//bRGjBhhdGpAk0GdV4k6rypqPXfUeu6o886POs8/6zyaazDU2LFjtW3bNm3atMnoVAx19dVXq7i4WFarVcuXL9eoUaP00UcfNcnCq6SkROPHj9f777+v4OBgo9PxCz//BaZHjx5KSkrSFVdcoddff13333+/gZkZx+FwqHfv3po2bZokqWfPntq2bZtycnIouiS9+OKLuuWWW9SuXTujUzHM66+/rn/+859asmSJrrnmGhUXF2vChAlq164d/4wADYQ6rxJ1njtqvaqo9dxR550fdZ5/1nk012CYcePGac2aNfr444/VoUMHo9MxVGBgoK688kpJUkJCgj799FPNnj1bL7zwgsGZNbzCwkIdPnxYvXr1co3Z7XZ9/PHHmjNnjioqKhQQEGBghsa77LLLdNVVV2nXrl1Gp2KYyMjIKv9R8otf/EJvvvmmQRn5j3379umDDz7QW2+9ZXQqhnr44Yc1adIk3X333ZKk7t27a9++fcrKyqIwBxoAdd451HnuqPUurKnXetR5NaPOq+SPdR7NNTQ4p9OpBx98UCtWrNCGDRua/OSU1XE4HKqoqDA6DUPcdNNN+vzzz93G0tPT1a1bN02cOLHJF1tS5QTAu3fv1r333mt0Koa59tprtWPHDrexr7/+WldccYVBGfmPl156SWFhYa4JXpuqkydPymx2n1o2ICBADofDoIyApoE678Kacp0nUevVRlOv9ajzakadV8kf6zyaa/Xshx9+cPvFYc+ePSouLlabNm3UsWNHAzMzztixY7VkyRK9/fbbatWqlUpLSyVJFotFLVq0MDi7hpeZmalbbrlFHTt21IkTJ7RkyRJt2LBB7733ntGpGaJVq1ZV5mW55JJL1LZt2yY7X8tf/vIXDR06VFdccYUOHjyoKVOmKCAgQMOHDzc6NcM89NBD6tevn6ZNm6a77rpLBQUFmj9/vubPn290aoZyOBx66aWXNGrUKDVr1rT/FT906FA9/fTT6tixo6655hpt3bpVs2bN0ujRo41ODRcZaj131HnuqPOqotarilrPHXVe9ajzzvHLOs+JerV+/XqnpCrbqFGjjE7NMNVdD0nOl156yejUDDF69GjnFVdc4QwMDHRefvnlzptuusm5bt06o9PyKwMGDHCOHz/e6DQMk5aW5oyMjHQGBgY627dv70xLS3Pu2rXL6LQMt3r1amdsbKwzKCjI2a1bN+f8+fONTslw7733nlOSc8eOHUanYjibzeYcP368s2PHjs7g4GBn586dnX/729+cFRUVRqeGiwy1njvqPHfUebVDrUet99+o86qizjvHH+s8k9PpdDZcKw8AAAAAAAC4eJgvHAIAAAAAAACgOjTXAAAAAAAAAA/RXAMAAAAAAAA8RHMNAAAAAAAA8BDNNQAAAAAAAMBDNNcAAAAAAAAAD9FcAwAAAAAAADxEcw0AAAAAAADwEM01AAAAAAAAwEM01wAAAAAAAAAP0VwDAAAAAAAAPERzDQAAAAAAAPAQzTUAOI+9e/fKZDJp4MCBstlsysjIUKdOndS8eXNNmDDB6PQAAADgIeo8AL7SzOgEAKAx+PHHHzVgwADt27dPAwYMUK9evdS6dWuj0wIAAICXqPMAeIvmGgDUQkFBgZKTk/XNN9/osssuMzodAAAA+Ah1HgBv8VgoANTS//7v/1JwAQAAXISo8wB4g+YaANRCZGSkevfubXQaAAAA8DHqPADeorkGALXQsWNHo1MAAABAPaDOA+AtmmsAUAvBwcFGpwAAAIB6QJ0HwFs01wAAAAAAAAAP0VwDAAAAAAAAPERzDQAAAAAAAPAQzTUAAAAAAADAQzTXAAAAAAAAAA+ZnE6n0+gkAAAAAAAAgMaIO9cAAAAAAAAAD9FcAwAAAAAAADxEcw0AAAAAAADwEM01AAAAAAAAwEM01wAAAAAAAAAP0VwDAAAAAAAAPERzDQAAAAAAAPAQzTUAAAAAAADAQzTXAAAAAAAAAA/RXAMAAAAAAAA8RHMNAAAAAAAA8BDNNQAAAAAAAMBD/x8yr9tA1/+46QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 1500x400 with 2 Axes>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "plt.figure(figsize=(15, 4))\n",
    "plt.subplot(1,2,1)\n",
    "plt.scatter(r, p)\n",
    "plt.xlabel(\"r\", fontsize=15)\n",
    "plt.ylabel(\"Probability\", fontsize=15)\n",
    "\n",
    "plt.subplot(1,2,2)\n",
    "plt.scatter(r, distances, color=\"#FF6600\")\n",
    "plt.xlabel(\"r\", fontsize=15)\n",
    "plt.ylabel(\"Distance\", fontsize=15)\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "id": "57898331",
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.13"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 5
}