BOCS_MAP.ipynb

bocs_map.ipynb

{
  "cells": [
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "view-in-github",
        "colab_type": "text"
      },
      "source": [
        "<a href=\"https://colab.research.google.com/gist/mohzeki222/9ca4cae7191de9dcc4da99d40804920e/bocs_map.ipynb\" target=\"_parent\"><img src=\"https://colab.research.google.com/assets/colab-badge.svg\" alt=\"Open In Colab\"/></a>"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "V1CMs3ILjIBo"
      },
      "source": [
        "量子アニーリングマシン（D-Waveマシン）は組合せ最適化問題を解くために用いられますが，そのためにはQUBOという問題形式に書き換える必要がありました．ですが，実際の世の中の問題全てがQUBOで表される訳ではありません．そのようなときに用いられるのがブラックボックス最適化と呼ばれる手法です．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "hBtTi45qNC7e"
      },
      "source": [
        "まずはいつものようにD-Waveマシンを使うためのライブラリをインストールしておきましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "fRMSRlutNCsX",
        "outputId": "c6c38755-e715-49d6-81ee-c8bd279ea44c"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting dwave-ocean-sdk\n",
            "  Downloading dwave_ocean_sdk-6.9.0-py3-none-any.whl (8.4 kB)\n",
            "Collecting dimod==0.12.14 (from dwave-ocean-sdk)\n",
            "  Downloading dimod-0.12.14-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (18.7 MB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m18.7/18.7 MB\u001b[0m \u001b[31m35.8 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting dwave-cloud-client==0.11.3 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_cloud_client-0.11.3-py3-none-any.whl (138 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m138.7/138.7 kB\u001b[0m \u001b[31m4.9 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting dwave-greedy==0.3.0 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_greedy-0.3.0-py3-none-any.whl (10 kB)\n",
            "Collecting dwave-hybrid==0.6.11 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_hybrid-0.6.11-py3-none-any.whl (77 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m77.1/77.1 kB\u001b[0m \u001b[31m6.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting dwave-inspector==0.4.4 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_inspector-0.4.4-py3-none-any.whl (31 kB)\n",
            "Collecting dwave-neal==0.6.0 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_neal-0.6.0-py3-none-any.whl (8.7 kB)\n",
            "Collecting dwave-networkx==0.8.14 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_networkx-0.8.14-py3-none-any.whl (102 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m102.4/102.4 kB\u001b[0m \u001b[31m11.4 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting dwave-preprocessing==0.6.5 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_preprocessing-0.6.5-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (3.4 MB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m3.4/3.4 MB\u001b[0m \u001b[31m23.2 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting dwave-samplers==1.2.0 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_samplers-1.2.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (6.7 MB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m6.7/6.7 MB\u001b[0m \u001b[31m51.9 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting dwave-system==1.23.0 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_system-1.23.0-py3-none-any.whl (103 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m103.2/103.2 kB\u001b[0m \u001b[31m11.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting dwave-tabu==0.5.0 (from dwave-ocean-sdk)\n",
            "  Downloading dwave_tabu-0.5.0-py3-none-any.whl (9.2 kB)\n",
            "Collecting dwavebinarycsp==0.3.0 (from dwave-ocean-sdk)\n",
            "  Downloading dwavebinarycsp-0.3.0-py3-none-any.whl (35 kB)\n",
            "Collecting minorminer==0.2.13 (from dwave-ocean-sdk)\n",
            "  Downloading minorminer-0.2.13-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (10.3 MB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m10.3/10.3 MB\u001b[0m \u001b[31m58.2 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hCollecting penaltymodel==1.1.0 (from dwave-ocean-sdk)\n",
            "  Downloading penaltymodel-1.1.0-py3-none-any.whl (36 kB)\n",
            "Requirement already satisfied: numpy<2.0.0,>=1.17.3 in /usr/local/lib/python3.10/dist-packages (from dimod==0.12.14->dwave-ocean-sdk) (1.25.2)\n",
            "Requirement already satisfied: requests[socks]>=2.18 in /usr/local/lib/python3.10/dist-packages (from dwave-cloud-client==0.11.3->dwave-ocean-sdk) (2.31.0)\n",
            "Requirement already satisfied: pydantic<3,>=2 in /usr/local/lib/python3.10/dist-packages (from dwave-cloud-client==0.11.3->dwave-ocean-sdk) (2.6.3)\n",
            "Collecting homebase>=1.0 (from dwave-cloud-client==0.11.3->dwave-ocean-sdk)\n",
            "  Downloading homebase-1.0.1-py2.py3-none-any.whl (11 kB)\n",
            "Requirement already satisfied: click>=7.0 in /usr/local/lib/python3.10/dist-packages (from dwave-cloud-client==0.11.3->dwave-ocean-sdk) (8.1.7)\n",
            "Requirement already satisfied: python-dateutil>=2.7 in /usr/local/lib/python3.10/dist-packages (from dwave-cloud-client==0.11.3->dwave-ocean-sdk) (2.8.2)\n",
            "Collecting plucky>=0.4.3 (from dwave-cloud-client==0.11.3->dwave-ocean-sdk)\n",
            "  Downloading plucky-0.4.3-py2.py3-none-any.whl (10 kB)\n",
            "Collecting diskcache>=5.2.1 (from dwave-cloud-client==0.11.3->dwave-ocean-sdk)\n",
            "  Downloading diskcache-5.6.3-py3-none-any.whl (45 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m45.5/45.5 kB\u001b[0m \u001b[31m5.3 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hRequirement already satisfied: packaging>=19 in /usr/local/lib/python3.10/dist-packages (from dwave-cloud-client==0.11.3->dwave-ocean-sdk) (23.2)\n",
            "Requirement already satisfied: werkzeug>=2.2 in /usr/local/lib/python3.10/dist-packages (from dwave-cloud-client==0.11.3->dwave-ocean-sdk) (3.0.1)\n",
            "Requirement already satisfied: typing-extensions>=4.5.0 in /usr/local/lib/python3.10/dist-packages (from dwave-cloud-client==0.11.3->dwave-ocean-sdk) (4.10.0)\n",
            "Collecting authlib<2,>=1.2 (from dwave-cloud-client==0.11.3->dwave-ocean-sdk)\n",
            "  Downloading Authlib-1.3.0-py2.py3-none-any.whl (223 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m223.7/223.7 kB\u001b[0m \u001b[31m17.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hRequirement already satisfied: importlib-metadata>=5.0.0 in /usr/local/lib/python3.10/dist-packages (from dwave-cloud-client==0.11.3->dwave-ocean-sdk) (7.0.1)\n",
            "Requirement already satisfied: networkx in /usr/local/lib/python3.10/dist-packages (from dwave-hybrid==0.6.11->dwave-ocean-sdk) (3.2.1)\n",
            "Requirement already satisfied: Flask>=2.2 in /usr/local/lib/python3.10/dist-packages (from dwave-inspector==0.4.4->dwave-ocean-sdk) (2.2.5)\n",
            "Requirement already satisfied: scipy>=1.7.3 in /usr/local/lib/python3.10/dist-packages (from dwave-system==1.23.0->dwave-ocean-sdk) (1.11.4)\n",
            "Collecting fasteners>=0.15 (from minorminer==0.2.13->dwave-ocean-sdk)\n",
            "  Downloading fasteners-0.19-py3-none-any.whl (18 kB)\n",
            "Collecting rectangle-packer>=2.0.1 (from minorminer==0.2.13->dwave-ocean-sdk)\n",
            "  Downloading rectangle_packer-2.0.2-cp310-cp310-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl (305 kB)\n",
            "\u001b[2K     \u001b[90m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m305.6/305.6 kB\u001b[0m \u001b[31m24.6 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n",
            "\u001b[?25hRequirement already satisfied: cryptography in /usr/local/lib/python3.10/dist-packages (from authlib<2,>=1.2->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (42.0.5)\n",
            "Requirement already satisfied: Jinja2>=3.0 in /usr/local/lib/python3.10/dist-packages (from Flask>=2.2->dwave-inspector==0.4.4->dwave-ocean-sdk) (3.1.3)\n",
            "Requirement already satisfied: itsdangerous>=2.0 in /usr/local/lib/python3.10/dist-packages (from Flask>=2.2->dwave-inspector==0.4.4->dwave-ocean-sdk) (2.1.2)\n",
            "Requirement already satisfied: zipp>=0.5 in /usr/local/lib/python3.10/dist-packages (from importlib-metadata>=5.0.0->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (3.17.0)\n",
            "Requirement already satisfied: annotated-types>=0.4.0 in /usr/local/lib/python3.10/dist-packages (from pydantic<3,>=2->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (0.6.0)\n",
            "Requirement already satisfied: pydantic-core==2.16.3 in /usr/local/lib/python3.10/dist-packages (from pydantic<3,>=2->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (2.16.3)\n",
            "Requirement already satisfied: six>=1.5 in /usr/local/lib/python3.10/dist-packages (from python-dateutil>=2.7->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (1.16.0)\n",
            "Requirement already satisfied: charset-normalizer<4,>=2 in /usr/local/lib/python3.10/dist-packages (from requests[socks]>=2.18->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (3.3.2)\n",
            "Requirement already satisfied: idna<4,>=2.5 in /usr/local/lib/python3.10/dist-packages (from requests[socks]>=2.18->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (3.6)\n",
            "Requirement already satisfied: urllib3<3,>=1.21.1 in /usr/local/lib/python3.10/dist-packages (from requests[socks]>=2.18->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (2.0.7)\n",
            "Requirement already satisfied: certifi>=2017.4.17 in /usr/local/lib/python3.10/dist-packages (from requests[socks]>=2.18->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (2024.2.2)\n",
            "Requirement already satisfied: PySocks!=1.5.7,>=1.5.6 in /usr/local/lib/python3.10/dist-packages (from requests[socks]>=2.18->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (1.7.1)\n",
            "Requirement already satisfied: MarkupSafe>=2.1.1 in /usr/local/lib/python3.10/dist-packages (from werkzeug>=2.2->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (2.1.5)\n",
            "Requirement already satisfied: cffi>=1.12 in /usr/local/lib/python3.10/dist-packages (from cryptography->authlib<2,>=1.2->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (1.16.0)\n",
            "Requirement already satisfied: pycparser in /usr/local/lib/python3.10/dist-packages (from cffi>=1.12->cryptography->authlib<2,>=1.2->dwave-cloud-client==0.11.3->dwave-ocean-sdk) (2.21)\n",
            "Installing collected packages: rectangle-packer, plucky, homebase, fasteners, diskcache, dimod, penaltymodel, dwave-samplers, dwave-preprocessing, dwave-networkx, minorminer, dwavebinarycsp, dwave-tabu, dwave-neal, dwave-greedy, authlib, dwave-cloud-client, dwave-system, dwave-inspector, dwave-hybrid, dwave-ocean-sdk\n",
            "Successfully installed authlib-1.3.0 dimod-0.12.14 diskcache-5.6.3 dwave-cloud-client-0.11.3 dwave-greedy-0.3.0 dwave-hybrid-0.6.11 dwave-inspector-0.4.4 dwave-neal-0.6.0 dwave-networkx-0.8.14 dwave-ocean-sdk-6.9.0 dwave-preprocessing-0.6.5 dwave-samplers-1.2.0 dwave-system-1.23.0 dwave-tabu-0.5.0 dwavebinarycsp-0.3.0 fasteners-0.19 homebase-1.0.1 minorminer-0.2.13 penaltymodel-1.1.0 plucky-0.4.3 rectangle-packer-2.0.2\n"
          ]
        }
      ],
      "source": [
        "!pip install dwave-ocean-sdk"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "k4oeMQB0bZvo"
      },
      "source": [
        "あとは今回は最適化の様子を可視化するために必要なモジュールをインストールしておきましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "d56TFLl8cCtm",
        "outputId": "cf3029fd-89cc-4e91-c238-f12d4f5ef075"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "Collecting APNG\n",
            "  Downloading apng-0.3.4-py2.py3-none-any.whl (8.2 kB)\n",
            "Installing collected packages: APNG\n",
            "Successfully installed APNG-0.3.4\n"
          ]
        }
      ],
      "source": [
        "!pip install APNG"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "FzpGy4ehfZS7"
      },
      "source": [
        "そして今回必要なライブラリも予め導入しておきましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "FUrDYtmFfW49"
      },
      "outputs": [],
      "source": [
        "import numpy as np\n",
        "import matplotlib.pyplot as plt\n",
        "from sklearn.linear_model import LinearRegression, Ridge\n",
        "from tqdm import tqdm\n",
        "from neal import SimulatedAnnealingSampler\n",
        "from IPython.display import display, clear_output, Image\n",
        "from apng import APNG"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "5Kgas93zhUzd"
      },
      "source": [
        "ブラックボックス最適化についてもう少し具体的に説明します．  \n",
        "問題設定として，データ$\\vec{x}=(x_1, x_2, ..., x_N)$が与えられたときに解$y=f(\\vec{x})$で出力されるとします．この関数$f$の中身がわからないときに関数$f$をブラックボックス関数といいます．  \n",
        "今回はランダムなQUBOをブラックボックス関数とします．ここで重要なのは**この目的関数は必ずしもQUBOの形である必要はありません**．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Qda9cbBooAIl"
      },
      "source": [
        "## ブラックボックス関数"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "GNMmBC_poGp7"
      },
      "source": [
        "ブラックボックス関数としては，ランダムなQUBOを用意しましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "uG6uj8vSpkXh"
      },
      "outputs": [],
      "source": [
        "N = 32\n",
        "QUBO = np.random.randn(N, N)\n",
        "QUBO = np.triu(QUBO)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "MalmB8ECqOVB"
      },
      "source": [
        "このQUBOは以下のような式で表されます．\n",
        "\\begin{equation}\n",
        "  {\\rm QUBO} = \\sum_{i=1}^N \\sum_{j=i}^N Q_{ij}x_i x_j, \\quad Q_{ij} \\sim N(0, 1)\n",
        "\\end{equation}\n",
        "$Q_{ij} \\sim N(0, 1)$は$Q_{ij}$が平均0，標準偏差1の標準正規分布に従う乱数であることを示します．  \n",
        "このQUBOが実際どのような形になっているのか見てみましょう．ついでに，このQUBOを可視化するというのはこの先何回も行うので関数で定義しておくと便利です．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "CSthTgVqqVr4"
      },
      "outputs": [],
      "source": [
        "def visualize_Q(Q, cmap=\"bwr\"):\n",
        "  cmap_range = np.abs(Q).max()\n",
        "  plt.imshow(Q, cmap=cmap, vmin=-cmap_range, vmax=cmap_range)\n",
        "  plt.colorbar()\n",
        "  plt.show()"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1TahZIaYu3oT"
      },
      "source": [
        "それでは実際に可視化してみましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 430
        },
        "id": "L4VdXalsqlxL",
        "outputId": "dc41a281-7640-45a0-8e1d-5c80eb181d71"
      },
      "outputs": [
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 2 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAeoAAAGdCAYAAADdSjBDAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAAAy80lEQVR4nO3deXhVVZb38V+CJIwJhsxNwADKUAzaUWIckCEy+BZCgTZWlSUgL5RUYhWGsiT9KjiUHad2bAqtbgZtxQFLROgWC0FC0w0okRTgkDY01QlCgkCTGwJJkJz3D5uUgQB7Jzfcc3K+n+c5z0NuVtZdx4tZrH3P3SfMcRxHAADAlcJDXQAAADg7GjUAAC5GowYAwMVo1AAAuBiNGgAAF6NRAwDgYjRqAABcjEYNAICLXRTqAk5XV1enffv2qXPnzgoLCwt1OQAAS47jqLKyUsnJyQoPb7l5sLq6WrW1tc3OExERoXbt2gWhopbhuka9b98+paSkhLoMAEAzlZaWqlu3bi2Su7q6Wqnt26ssCLkSExO1Z88e1zbrFmvUCxYs0JNPPqmysjINHjxYL7zwgoYMGXLen+vcubMkqXToUEVdZFjec8+ZF7Zpk3mspIOjfmIcG7t8oVXur2+eZRz78stWqTVsmHnsNRl2u8iu/dBupePIEfPYLl2sUuvGXv9lHvzZZ3bJMzPNY/fts0q96vNeVvHjMg6aB0dHW+WurG5rHPvnP1ulVkGBeezU/h/bJf/0U/PYiy+2yz1qlFV4TQfz/JGP3G9XS06Oceimz2OsUm/fbh77f/6PeezRowFdf31K/e/zllBbW6sySaVhYYpqRp6ApJSyMtXW1vqrUb/55pvKycnRiy++qPT0dD377LMaPXq0ioqKFB8ff86fPbXcHXXRReaN2uYvQ/v25rGSajqb/xWIsnyRAxa5IyOtUqtjR/PYqCi7Rt2hg12jtlmZ6tDBKrWibF576+QW//sHAlapO3Sw+9US1bnGItgud1iEeaPu1MkqtdX/blE2f2ltk1v+f2/737DG4vWMtP2f2eLveMeOdnXb/MpqSs+9EG9fRkmKas7zeOB2Fy3y5sHTTz+tGTNmaNq0aerfv79efPFFdejQQYsXL26JpwMA+FV4ePMPCwsXLtSgQYMUFRWlqKgoZWRk6P3332+hk/tO0Bt1bW2tCgoKlPm9ZcPw8HBlZmZq8+bNZ8TX1NQoEAg0OAAAMHKBG3W3bt302GOPqaCgQNu2bdOIESM0fvx4fWb71pqFoDfqgwcP6uTJk0pISGjweEJCgsrKznzbPy8vT9HR0fUHF5IBAIxd4EY9btw43XTTTbr00kt12WWX6dFHH1WnTp20ZcuWFjpBF3yOOjc3VxUVFfVHaWlpqEsCAPjM6Su7NTXnvy7k5MmTeuONN1RVVaWMjIwWqy3oF5PFxsaqTZs2Ki8vb/B4eXm5EhMTz4iPjIy0v7gCAADpu4k4CBeTnb6aO3/+fD344ION/sjOnTuVkZGh6upqderUSStWrFD//v2bXsN5BL1RR0REKC0tTevWrdOECRMkfbeJybp165SdnR3spwMA+FmQGnVpaamivne1/7kGyD59+qiwsFAVFRV6++23NWXKFOXn57dYs26Rj2fl5ORoypQpuvLKKzVkyBA9++yzqqqq0rRp01ri6QAAaJZTV3GbiIiIUO/evSVJaWlp+uSTT/Tcc8/ppZdeapHaWqRRT548Wd98843mzZunsrIyXX755VqzZs0ZF5gBANAsQZqom6Ours7oPe2marGdybKzs5u31P3OO8abDjgyf5HCIiKsyogLP2Qc+3KXX1nlnmCxgcCAAVapZbVr39tvW+Uee9NNVvEn25lvZHH4sFVqKedB89jHHrPLvWOHcejODulWqX9080m7WjZ9bhx658s3WKVefMULxrGDb7/dKnevXhY7gv1zoVVuzZxpHPrRf9hdB3O15QZVETa/SceMsUv+3nvGoUdi7FYtbX4d9qr7yjg24By1qqNZLnCjzs3N1dixY9W9e3dVVlZq2bJl2rBhgz744IOm13AertvrGwAAtzpw4IDuuOMO7d+/X9HR0Ro0aJA++OAD3XjjjS32nDRqAIB3XeCJetGiRU1/riaiUQMAvCsszHrTkgbq6oJXSwsJ+YYnAADg7JioAQDe1YRtQL2GRg0A8C4aNQAALuaDRt26zw4AAI9jogYAeJcPJmoaNQDAu2jU3hAmiw+s2+7F+eijxqFTxo2zy73+iHHorVEdrFKvKBxtHJtquaXhj2433xJUklbM224cG2eVWTqx+J+NY9tWV1rlnvai+bagd91llVon6tpYxbeNjzeOXfxMhV0xgQnmsV9/bZW6U8K3xrEn/u8sq9xtP91qHDt8+ctWuXX139vFf15sHnvJJXa5L7/cOLT0dbvUdx973Dh2b/v7jGMrvw3YFYJzahWNGgDgU0zUAAC4mA8ades+OwAAPI6JGgDgXT6YqGnUAADvau5NOSzvnhUKrfufIQAAeBwTNQDAu5q79O2BiZpGDQDwLho1AAAu5oNGzXvUAAC4GBM1AMC7fDBR+69RW+71PeXQ08axv4yyKyXtSou9pzt1sspd8rxF8C23WOVeYbmn+brD2caxI2P/ZJW77eFy8+ANG6xyz5s32Th21y6r1Jo40S5+VTeLF/TBB61y//blFOPY+7ustMr9xUjz1z7WKrMUd9ll5sGz7PYR16ZNVuE7E280jh2o3Xa1HDtmHJo9w+73hI5MNQ7tlnuncWygttaujubwQaNm6RsAABfz30QNAGg9fDBR06gBAN7lg0bN0jcAAC7GRA0A8C4fTNQ0agCAdzX3phx1dcGrpYWw9A0AgIsxUQMAvKu5S98euJc1jRoA4F00agAAXIxGjZeXWlwRuHevXfLDFhcx/Ou/WqX+VfWfzYP/5m+scuf3nm4Vf+SgRXDdAavcWrXKPDYhwSp1agfz7Um3Vdvlvv9+q3Dtjl1oHBt5wi5327YWwaNHW+XuF19hHlxcbJV7e3iacewVtf9plVvXXWcV/vl75rF/924vq9yv/90e49hZsyOtci98ynzL0VdHLDaOPX48IL32mlUtODsaNQDAu5ioAQBwMR80avdXCACAjzFRAwC8ywcTNY0aAOBdPmjU7q8QAAAfY6IGAHiXDyZqGjUAwLuae1OOsLDg1dJC3P9PCQAAfIyJGgDgXSx9AwDgYjRqWOnWzSr8t4+avzcyZ04Pq9xXX20e+6fsf7TKfcOB5VbxOnbMOPT4TVOsUi8tvtE4dtaHk61y2+z3fGui5V7SHS+2Cn9u9QDj2G+/tStlwgTz2EUbL7XKPX3YbuPYLzqY790tSVck/o958NFkq9wqK7MKn3xLd+PYmJg2LVbLrl2pVqn/+Z2OxrE/u+W4cWwgcFwzZ1qV0nQ+aNTurxAAAB8LeqN+8MEHFRYW1uDo27dvsJ8GAIC/TNTNOVyuRZa+f/CDH+jDDz/8y5NcxAo7AKAF+GDpu0U66EUXXaTExMSWSA0AgK+0yD8lvvrqKyUnJ6tnz5766U9/qpKSkrPG1tTUKBAINDgAADBygZe+8/LydNVVV6lz586Kj4/XhAkTVFRU1EIn952gN+r09HQtXbpUa9as0cKFC7Vnzx5df/31qqysbDQ+Ly9P0dHR9UdKSkqwSwIAtFYXuFHn5+crKytLW7Zs0dq1a3XixAmNGjVKVVVVLXSCLbD0PXbs2Po/Dxo0SOnp6erRo4feeustTZ8+/Yz43Nxc5eTk1H8dCARo1gAAV1qzZk2Dr5cuXar4+HgVFBRo6NChLfKcLX6VV5cuXXTZZZepuLi40e9HRkYqMjKypcsAALRGQbqY7PS3XU17U0VFhSQpJiam6TWcR4tf7nb06FHt3r1bSUlJLf1UAAC/OXVTjqYe/3tTjpSUlAZvw+bl5Z33qevq6jR79mxde+21GjDAfFMiW0GfqH/9619r3Lhx6tGjh/bt26f58+erTZs2+vGPfxzspwIAIChKS0sVFRVV/7XJNJ2VlaVdu3Zp06ZNLVla8Bv13r179eMf/1iHDh1SXFycrrvuOm3ZskVxcXHBfirPu///Ocaxzz1vdyu2P/1ykXnslTOscg8eZF63JGn9euPQ9p8XWKWe9cN48+Daa6xyr93bzzj2xi6fWOXWH/5gFf6rm8y3b3zgvauscveJ2m8c2/MOu5Wxf9vSyzj2+svKrXKv3JhgHDv+W/O/g5L0D/snWcVnd1xiHJs5dZpVbv2/Vcah//a3R6xSfxA+9vxBp1RXt0xscwVp6TsqKqpBoz6f7OxsrV69Whs3blQ3y+2jbQW9Ub/xxhvBTgkAQOMu8IYnjuPo7rvv1ooVK7Rhwwalptrtr94UbBkGAPCuC9yos7KytGzZMq1cuVKdO3dW2f/eNCU6Olrt27dveh3n4P690wAAcImFCxeqoqJCw4YNU1JSUv3x5ptvtthzMlEDALwrBEvfFxqNGgDgXT64KYf7KwQAwMeYqAEA3uWDiZpGDQDwLh80avdXCACAjzFRAwC8ywcTNY0aAOBdp27K0ZyfdzkatUf86peWn93bO8o4tOxzu9SBgN1f7KeeH2kca7sDbXuZ74GtH/7QKveNyx4xjl13zQNWuUfOtNwbuLDQODQ93S71/BfN9+9+6I7dVrmvvtp8r+9/22K+d7ckjf/WYr/0vXutck+w2+pbC1eZ79896+OtVrn/NPnvjGNPu1PjeY0+vNI49lDdeOPYSqeNXSE4Jxo1AMC7WPoGAMDFaNQAALiYDxq1+ysEAMDHmKgBAN7lg4maRg0A8C4fNGr3VwgAgI8xUQMAvMsHEzWNGgDgXT5o1O6vEAAAH2Oibq26mW9ROXqf3ZaGmj3bKnzAv242jt22za6U6zt9aR68eLFV7k/ueME4duTlJ6xya1+tXfyQIcahvQ7Ypf7hJ/ONY1/d/JBV7tt/Mdo49vrf/94q94mrzff5fOklq9RKsvxfYtbEcvPgfRFWuQcNMo8N2/RvVrm/6m++Lej6t83zHrfY2bfZfDBR06gBAN7lg5tyuP+fEgAA+BgTNQDAu1j6BgDAxWjUAAC4mA8atfsrBADAx5ioAQDe5YOJmkYNAPAuHzRq91cIAICPMVEDALzLBxM1jRoA4F00aviCxT7SkvTFYvO9uyUp3GLv6esvs9gzWZI69DaPff55q9RXFX9lHjxqplVuvW2xcbKk/wnvahzb71iBVe78Eeb7d9/eyS63pk41Dt28r4dV6ox25n9Xsqd1ssq9aXtHq3irX/YxMVapw45VmQdb7PEvSZf+1zrj2P9oN9I41nGsysB50KgBAN7FRA0AgItxUw4AABBKTNQAAO9i6RsAABejUQMA4GI+aNTurxAAAB9jogYAeJcPJmoaNQDAu3zQqN1fIQAAPsZEDQDwLh9M1DRqWOvX13Ij3zvvNA7Nn7rEKvUN15wwD66rs8p9suelxrFtLPdL319rvne3JCWVbTcPjoy0yt3uW/PYnRFpdrmvNI/P+HylVe47/3G8cezdd1ulVlKSXfxXR+KMYy/ttN8u+bZt5rF//rNV6u2DphjHXh5rnvfoUasymscHjdr9FQIA4GPWjXrjxo0aN26ckpOTFRYWpnfffbfB9x3H0bx585SUlKT27dsrMzNTX31lcRciAABMnZqom3O4nHWFVVVVGjx4sBYsWNDo95944gk9//zzevHFF7V161Z17NhRo0ePVnV1dbOLBQCggVM35Wjq4YGbcli/Rz127FiNHTu20e85jqNnn31W999/v8aP/+79o1deeUUJCQl69913ddtttzWvWgAAfCaoM/+ePXtUVlamzMzM+seio6OVnp6uzZs3N/ozNTU1CgQCDQ4AAIyw9G2nrKxMkpSQkNDg8YSEhPrvnS4vL0/R0dH1R0pKSjBLAgC0ZiFo1Oe7VivYQv5PidzcXFVUVNQfpaWloS4JAOAVIWjU57tWK9iC+jnqxMRESVJ5ebmSvvdBxPLycl1++eWN/kxkZKQiLT/3CQBAqJzrWq2WENSJOjU1VYmJiVq3bl39Y4FAQFu3blVGRkYwnwoAgKBN1KdfK1VTUxPiE/sL60Z99OhRFRYWqrCwUNJ3F5AVFhaqpKREYWFhmj17tn7729/qvffe086dO3XHHXcoOTlZEyZMCHLpAADfC1KjTklJaXC9VF5eXohP7C+sl763bdum4cOH13+dk5MjSZoyZYqWLl2q3/zmN6qqqtLMmTN15MgRXXfddVqzZo3atWsXvKrhLYsXG4e+/Uu71Dd0LTIP/ugjq9zvp5rvOxn9w8etcl8fa7H1qSQ1fi1m41JTrVLHHzCP3bTJKrV+NuJr49hvrjHfElSSFh8232725KBpVrkPHrQKt/rvcmnJG3bJk5PNY0eMsEp9xeKHzIPLy41DA7W1VnW4QWlpqaKiouq/dtNbstaNetiwYXKcs+/1HBYWpocfflgPP/xwswoDAOC8grTXd1RUVING7SbclAMA4F0+uCkHjRoAAAtHjx5VcXFx/denrtWKiYlR9+7dg/58NGoAgHeFYKI+37VawUajBgB4Vwga9fmu1Qo29y/OAwDgY0zUAADvOnWby+b8vMvRqAEA3sVV3wAAuJgPGrX7KwQAwMeYqAEA3uWDiZpGDVd54XnLjzz8doVx6OorHrBKvX69eezT/2O3l7S6zbMKf2T1FcaxDxy9zyp3+C/M9ym/5Rar1NJNtxuHljxltxd7nMWNftos/Aer3Lv6ZVvFT5pw0jj2pO6xyv3ee+axP9pm/v+DJGnDBuPQTY+Yvz5VVQFp0SK7WprKB43a/RUCAOBjTNQAAO/ywURNowYAeJcPGrX7KwQAwMeYqAEA3uWDiZpGDQDwLho1AAAu5oO9vt3/TwkAAHyMiRoA4F0sfQMA4GI0asDl7r/fOLS9xZagkvT0MxbvXb3yil1ySzffbBF85Car3AcOmMf2qNxllVs5OcahacVv2uVOHmoc+moXuy1B2x60K8UJb2Mc2+YfXrDKPSH7buPYf371R1a5f1b8S+PY666qMY4NBMxjcX40agCAdzFRAwDgYj5o1O6vEAAAH2OiBgB4lw8maho1AMC7fNCo3V8hAAA+xkQNAPAuH0zUNGoAgHfRqAEAcDFuygEAAEKJiRoA4F0sfQOtx8gRjt0PzJtnHFo54WdWqa+5xq6Ud981j/2TbrDKfVW3Q+bBJSescqu21jw2IsIq9cqPk4xjJ0ywSq1ORQVW8bv/K804tsMt5nt3S1LhGvPYnx18xiq3/vhH49BNn0Qax1ZVmcc2mw8atfsrBADAx5ioAQDe5YOJmkYNAPAuHzRq91cIAICPMVEDALzLBxM1jRoA4F0+aNTurxAAAB9jogYAeJcPJmoaNQDAu2jUAAC4mA9uykGjBs7mwQeNQ+sCdql3znzBKn57wHzbySu03Sr38Q5XGMe2b7vfKvc/Hp5kHDtjqt32pOO3bjIP/uavrHJb7dkq6dGvzbcQXdzvSavcr+pe49ixB/ZZ5da33xqHXpe02zg2UFlpVwfOiUYNAPAulr4BAHAxHzRq6wo3btyocePGKTk5WWFhYXr3tCWiqVOnKiwsrMExZsyYYNULAICvWE/UVVVVGjx4sO68805NnDix0ZgxY8ZoyZIl9V9HRl7AW54BAPzDBxO1daMeO3asxo4de86YyMhIJSYmNrkoAACM+KBRt0iFGzZsUHx8vPr06aNZs2bp0KGz35i+pqZGgUCgwQEAAL4T9EY9ZswYvfLKK1q3bp0ef/xx5efna+zYsTp58mSj8Xl5eYqOjq4/UlJSgl0SAKC1OjVRN+dwuaBXeNttt+nmm2/WwIEDNWHCBK1evVqffPKJNmzY0Gh8bm6uKioq6o/S0tJglwQAaK1C1KgXLFigSy65RO3atVN6ero+/vjjIJ/YX7T4PyV69uyp2NhYFRcXN/r9yMhIRUVFNTgAAHCrN998Uzk5OZo/f74+/fRTDR48WKNHj9aBAwda5PlavFHv3btXhw4dUlJSUks/FQDAb0IwUT/99NOaMWOGpk2bpv79++vFF19Uhw4dtHjx4hY4wSZc9X306NEG0/GePXtUWFiomJgYxcTE6KGHHtKkSZOUmJio3bt36ze/+Y169+6t0aNHB7VwAACCddX36RcyR0ZGNvrR4traWhUUFCg3N/d7KcKVmZmpzZs3N72Oc7Bu1Nu2bdPw4cPrv87JyZEkTZkyRQsXLtSOHTv08ssv68iRI0pOTtaoUaP0yCOP8FlqtGrRUY5V/MlfmO/dLUlXvP6qcezOwbdb5R549WDz4FfN65Cks2y10Ki95W2tcv/uX68zjv27Rxq/mPVsTj74iFX8iNfNY//Q3nzvbkm63OLdwIIYu33EAwfNY+vqzGOrqi7cp3cchclR02+scepnT7+Qef78+Xqwkf3+Dx48qJMnTyohIaHB4wkJCfryyy+bXMe5WDfqYcOGyXHO/kvpgw8+aFZBAABcaKWlpQ2ukXLTcMle3wAAz6qrs5v2G/t5ScYXM8fGxqpNmzYqLy9v8Hh5eXmLbfTl/g+QAQBwFqcadXMOGxEREUpLS9O6deu+V0Od1q1bp4yMjCCf3XeYqAEAsJCTk6MpU6boyiuv1JAhQ/Tss8+qqqpK06ZNa5Hno1EDADwrWEvfNiZPnqxvvvlG8+bNU1lZmS6//HKtWbPmjAvMgoVGDQDwrFA0aknKzs5WdnZ205/YAu9RAwDgYkzUAADPCtVEfSHRqAEAnuWHRs3SNwAALsZEDQDwLD9M1DRqIATahNvtDS6LrXnrBtrt9a3HHjMOXflfA61Sf/q2eWz//lapddNNFsF5eVa53+13v1X87Z1XGsf+4dvxVrlv/PIF49gnq+32kLe5aLm9jhvHBgLmsc1FowYAwMUcp3nN9hy3rnAN3qMGAMDFmKgBAJ7F0jcAAC7mh0bN0jcAAC7GRA0A8Cw/TNQ0agCAZ/mhUbP0DQCAizFRAwA8yw8TNY0aAOBZNGoA7vDKK8ahhzdY5h4xwjh0/H/utEo9fkiseXC43TtxlR0SjGNf/bPllqCX77aKX/2F+bagxV9apdYnI8y3Bb33ou12ydXXvI5d7Y1jjx49YVkHzoVGDQDwLCZqAABcjEYNAICLcVMOAAAQUkzUAADPYukbAAAX80OjZukbAAAXY6IGAHiWHyZqGjUAwLP80KhZ+gYAwMWYqAEAnuWHiZpGDbQyw4dZ7uDwcaFx6J74dKvUqbfdYBxb/la+Ve5PN5nHXnqpVWrp3nutwn/42mvGsQMHmu+ZLUnffmse++Q7V1jlvrdnhXHsVZccNY4NVFZa1dEcfmjULH0DAOBiTNQAAM/yw0RNowYAeBaNGgAAF+OmHAAAIKSYqAEAnsXSNwAALuaHRs3SNwAALsZEDQDwLD9M1DRqAIBn0agBtH5DhhiHppb8t13ud94xDk14+yWr1JcM/blxbL+sEVa5FR1tF19SYh7bro9V6l6rnzOOvbdfL6vcmv0H49C1P1liHFtVFWlXB86JRg0A8Cw/TNRWF5Pl5eXpqquuUufOnRUfH68JEyaoqKioQUx1dbWysrLUtWtXderUSZMmTVJ5eXlQiwYAQPpLo27O4XZWjTo/P19ZWVnasmWL1q5dqxMnTmjUqFGqqqqqj7nnnnu0atUqLV++XPn5+dq3b58mTpwY9MIBAPADq6XvNWvWNPh66dKlio+PV0FBgYYOHaqKigotWrRIy5Yt04gR370ntGTJEvXr109btmzR1VdfHbzKAQC+x9L3eVRUfHcv05iYGElSQUGBTpw4oczMzPqYvn37qnv37tq8eXOjOWpqahQIBBocAACYYOn7HOrq6jR79mxde+21GjBggCSprKxMERER6tKlS4PYhIQElZWVNZonLy9P0dHR9UdKSkpTSwIA+Mypm3I09WjVN+XIysrSrl279MYbbzSrgNzcXFVUVNQfpaWlzcoHAEBr0qSPZ2VnZ2v16tXauHGjunXrVv94YmKiamtrdeTIkQZTdXl5uRITExvNFRkZqchIPnMHALDHe9SncRxH2dnZWrFihdavX6/U1NQG309LS1Pbtm21bt26+seKiopUUlKijIyM4FQMAMD/8sN71FYTdVZWlpYtW6aVK1eqc+fO9e87R0dHq3379oqOjtb06dOVk5OjmJgYRUVF6e6771ZGRgZXfAMAfOXRRx/Vv/zLv6iwsFARERE6cuRIk/JYNeqFCxdKkoYNG9bg8SVLlmjq1KmSpGeeeUbh4eGaNGmSampqNHr0aP3ud79rUnEAAJyLm5e+a2trdeuttyojI0OLFi1qcp4wx3HXNW+BQEDR0dGqOHJEUVFRoS4HQDMcOhxmHNv1YNH5g76n4Kj5ntlpETutch9NHWgV3+kD8z2zt/ecZJW7ttY8NjnZKrXO8mGcRsXGmsdWVgY0eHC0KioqWuz3+Kle8eqrFerQoenPcexYQLff3rK1Ll26VLNnz74wEzUAAK3R6Xt4uOlC52ZteAIAQCgF62KylJSUBnt65OXlhfbEvoeJGgDgWcF6j7q0tLTB0vfZpum5c+fq8ccfP2fOL774Qn379m16UaehUQMAfC8qKsroPeo5c+bUXzx9Nj179gxSVd+hUQMAPOtCX/UdFxenuLi4pj9hE9CoAQCe5eaPZ5WUlOjw4cMqKSnRyZMnVVhYKEnq3bu3OnXqZJyHRg0A8Cw3N+p58+bp5Zdfrv/6iiuukCR99NFHZ+xHci5c9Q0AQAtYunSpHMc547Bp0hITNQDAw07d5rI5P+92NGoAgGe5eek7WGjUAFpM1xiLceXTEqvcaeF7jWO3x4y0yn1F+5NW8UdHm28L2uFrq9T63h2Dz+vAAbvcn35qHvvzyz4yjg1UVdkVgnOiUQMAPIuJGgAAF/NDo+aqbwAAXIyJGgDgWX6YqGnUAADP8kOjZukbAAAXY6IGAHiWHyZqGjUAwLNo1AAAuJgfGjXvUQMA4GJM1AAAz+KmHABwoWRmWoXv+izMOPaKf/yVXS1XX20V3ik+3ji2z4ABdrVYbPYdHh5plbpvX4vg2m/NY7+1iG0mlr4BAEBIMVEDADzLDxM1jRoA4Fl+aNQsfQMA4GJM1AAAz/LDRE2jBgB4lh8aNUvfAAC4GBM1AMCz/DBR06gBAJ5FowYAwMVo1ADgUgN+YL5J85LLn7PK/df97WoZ3LvKOHb7f3a0yt3usHlsvy9XWOXuUVJiHvyLX5jHBgJWdeDcaNQAAM/iphwAALiYH5a++XgWAAAuxkQNAPAsP0zUNGoAgGf5oVGz9A0AgIsxUQMAPMsPEzWNGgDgWX5o1Cx9AwDgYkzUAADP8sNETaMGAHgWjRoAWoFpUy33iXz2WavwJZ/eYxybnGxXyuhdf28cu7z7HKvcGZPMY7ttXGceXGW+93lz+aFR8x41AAAuZtWo8/LydNVVV6lz586Kj4/XhAkTVFRU1CBm2LBhCgsLa3DcddddQS0aAADpLzflaOrhhZtyWDXq/Px8ZWVlacuWLVq7dq1OnDihUaNGqeq0ZY4ZM2Zo//799ccTTzwR1KIBAJCa16Sbu2x+oVi9R71mzZoGXy9dulTx8fEqKCjQ0KFD6x/v0KGDEhMTg1MhAAA+1qz3qCsqKiRJMTExDR5/7bXXFBsbqwEDBig3N1fHjh07a46amhoFAoEGBwAAJpioz6Gurk6zZ8/WtddeqwEDBtQ//pOf/EQ9evRQcnKyduzYofvuu09FRUV65513Gs2Tl5enhx56qKllAAB8zA9XfTe5UWdlZWnXrl3atGlTg8dnzpxZ/+eBAwcqKSlJI0eO1O7du9WrV68z8uTm5ionJ6f+60AgoJSUlKaWBQBAq9KkRp2dna3Vq1dr48aN6tat2zlj09PTJUnFxcWNNurIyEhFRkY2pQwAgM8xUZ/GcRzdfffdWrFihTZs2KDU1NTz/kxhYaEkKSkpqUkFAgBwNjTq02RlZWnZsmVauXKlOnfurLKyMklSdHS02rdvr927d2vZsmW66aab1LVrV+3YsUP33HOPhg4dqkGDBrXICQAA0JpZXfW9cOFCVVRUaNiwYUpKSqo/3nzzTUlSRESEPvzwQ40aNUp9+/bVnDlzNGnSJK1atapFigcA+Jtbr/r+85//rOnTpys1NVXt27dXr169NH/+fNXW1lrnsl76PpeUlBTl5+dbFwEArjJ7tlX4tCP/Yxy758jFVrn3DzLfv/vWgtVWubU/wTz2e5/uOa/KSrs6msGtS99ffvml6urq9NJLL6l3797atWuXZsyYoaqqKj311FNWubgpBwDAs9zaqMeMGaMxY8bUf92zZ08VFRVp4cKFNGoAAGydvtlWS3wiqaKi4owNwkxw9ywAgGcF66YcKSkpio6Orj/y8vKCWmdxcbFeeOEF/fznP7f+WSZqAIBnBWvpu7S0VFFRUfWPn22anjt3rh5//PFz5vziiy/Ut2/f+q+//vprjRkzRrfeeqtmzJhhXSONGgDge1FRUQ0a9dnMmTNHU6dOPWdMz5496/+8b98+DR8+XNdcc41+//vfN6k2GjUAwLMu9MVkcXFxiouLM4r9+uuvNXz4cKWlpWnJkiUKD2/au800agCAZ7n1qu+vv/5aw4YNU48ePfTUU0/pm2++qf+e7W2gadQAAATZ2rVrVVxcrOLi4jPuiXG+PUlOx1XfAADPcuvOZFOnTpXjOI0etpioAQCe5dal72CiUQNAc3XpYhxad9gudVL1HuPYXZf80Cr3gLhy8+C1a81jjx+3qgPnRqMGAHgWEzUAAC5GowYAwMX80Ki56hsAABdjogYAeNapm3I05+fdjkYNAPAslr4BAEBIMVEDADzLDxM1jRoA4Fl+aNQsfQMA4GJM1AAAz/LDRE2jBoALqFdPy88DvfWxceilE1Ltcv/H5+axP/2peWwgIM2caVdLE/mhUbP0DQCAizFRAwA8yw8TNY0aAOBZNGoAAFzMD42a96gBAHAxJmoAgGdxUw4AAFysrk4KC2vez7sdS98AALgYEzUAwLP8MFHTqAEAnkWjBgCE1t/8jXHosqV2qQcMGG4cm7jXPG9lZTM6J85AowYAeBYTNQAALuaHRs1V3wAAuBgTNQDAs/wwUdOoAQCeRaMGAMDF/NCoeY8aAAAXY6IGAHgWN+UAAMDFmrt0zdI3AABoFiZqAIBn+WGiplEDQCsxbardG66bt5hfLn3smHne48etymgWPzRqlr4BAHAxq0a9cOFCDRo0SFFRUYqKilJGRobef//9+u9XV1crKytLXbt2VadOnTRp0iSVl5cHvWgAAKTvJuLmHm5n1ai7deumxx57TAUFBdq2bZtGjBih8ePH67PPPpMk3XPPPVq1apWWL1+u/Px87du3TxMnTmyRwgEA8EOjtnqPety4cQ2+fvTRR7Vw4UJt2bJF3bp106JFi7Rs2TKNGDFCkrRkyRL169dPW7Zs0dVXXx28qgEA8IkmX0x28uRJLV++XFVVVcrIyFBBQYFOnDihzMzM+pi+ffuqe/fu2rx581kbdU1NjWpqauq/DgQCTS0JAOAzXEzWiJ07d6pTp06KjIzUXXfdpRUrVqh///4qKytTRESEunTp0iA+ISFBZWVlZ82Xl5en6Ojo+iMlJcX6JAAA/uSHpW/rRt2nTx8VFhZq69atmjVrlqZMmaLPP/+8yQXk5uaqoqKi/igtLW1yLgCAv/ihUVsvfUdERKh3796SpLS0NH3yySd67rnnNHnyZNXW1urIkSMNpury8nIlJiaeNV9kZKQiIyPtKwcAwAea/Tnquro61dTUKC0tTW3bttW6devqv1dUVKSSkhJlZGQ092kAADjDqZtyNPVodTflyM3N1dixY9W9e3dVVlZq2bJl2rBhgz744ANFR0dr+vTpysnJUUxMjKKionT33XcrIyODK74BAC2iufejbnWN+sCBA7rjjju0f/9+RUdHa9CgQfrggw904403SpKeeeYZhYeHa9KkSaqpqdHo0aP1u9/9rkUKBwA0T8bV5l1qyVLzbnghtxD1A6tGvWjRonN+v127dlqwYIEWLFjQrKIAADDh5on65ptvVmFhoQ4cOKCLL75YmZmZevzxx5WcnGyVh72+AQCe5earvocPH6633npLRUVF+sMf/qDdu3frlltusc7D3bMAAGgB99xzT/2fe/Tooblz52rChAk6ceKE2rZta5yHRg0A8KxgLX2fvitmsD86fPjwYb322mu65pprrJq0xNI3AMDDgrX0nZKS0mCXzLy8vKDUd99996ljx47q2rWrSkpKtHLlSuscNGoAgO+VlpY22CUzNze30bi5c+cqLCzsnMeXX35ZH3/vvfdq+/bt+uMf/6g2bdrojjvukGN5BRtL3wAAzwrW0ndUVJSioqLOGz9nzhxNnTr1nDE9e/as/3NsbKxiY2N12WWXqV+/fkpJSdGWLVusNgKjUQMAPOtCfzwrLi5OcXFxTXquuv9dZ//+HSNN0KgBAJ7l1s9Rb926VZ988omuu+46XXzxxdq9e7ceeOAB9erVy3pbbd6jBgAgyDp06KB33nlHI0eOVJ8+fTR9+nQNGjRI+fn51leTu26iPvUm++mXygMAQuf4cZstRL/7/W170VRTOE6gmVNxy/SagQMHav369UHJ5bpGXVlZKUlK6d49xJUAAJqjsrJS0dHRLZI7IiJCiYmJKitLaXauxMRERUREBKGqlhHmXIh/8lioq6vTvn371LlzZ4V9742HQCCglJQUlZaWGl2Z51WcZ+vhh3OUOM/WJhjn6TiOKisrlZycrPDwlnuHtbq6WrW1tc3OExERoXbt2gWhopbhuok6PDxc3bp1O+v3TS+h9zrOs/XwwzlKnGdr09zzbKlJ+vvatWvn6gYbLFxMBgCAi9GoAQBwMc806sjISM2fPz+om6S7EefZevjhHCXOs7Xxy3l6iesuJgMAAH/hmYkaAAA/olEDAOBiNGoAAFyMRg0AgIt5plEvWLBAl1xyidq1a6f09HR9/PHHoS4pqB588MEzbj7et2/fUJfVLBs3btS4ceOUnJyssLAwvfvuuw2+7ziO5s2bp6SkJLVv316ZmZn66quvQlNsM5zvPKdOnXrGaztmzJjQFNtEeXl5uuqqq9S5c2fFx8drwoQJKioqahBTXV2trKwsde3aVZ06ddKkSZNUXl4eooqbxuQ8hw0bdsbredddd4Wo4qZZuHChBg0aVL+pSUZGht5///3677eG17I18USjfvPNN5WTk6P58+fr008/1eDBgzV69GgdOHAg1KUF1Q9+8APt37+//ti0aVOoS2qWqqoqDR48WAsWLGj0+0888YSef/55vfjii9q6das6duyo0aNHq7q6+gJX2jznO09JGjNmTIPX9vXXX7+AFTZffn6+srKytGXLFq1du1YnTpzQqFGjVFVVVR9zzz33aNWqVVq+fLny8/O1b98+TZw4MYRV2zM5T0maMWNGg9fziSeeCFHFTdOtWzc99thjKigo0LZt2zRixAiNHz9en332maTW8Vq2Ko4HDBkyxMnKyqr/+uTJk05ycrKTl5cXwqqCa/78+c7gwYNDXUaLkeSsWLGi/uu6ujonMTHRefLJJ+sfO3LkiBMZGem8/vrrIagwOE4/T8dxnClTpjjjx48PST0t5cCBA44kJz8/33Gc7167tm3bOsuXL6+P+eKLLxxJzubNm0NVZrOdfp6O4zg33HCD86tf/Sp0RbWQiy++2Pmnf/qnVvtaepnrJ+ra2loVFBQoMzOz/rHw8HBlZmZq8+bNIaws+L766islJyerZ8+e+ulPf6qSkpJQl9Ri9uzZo7Kysgava3R0tNLT01vd6ypJGzZsUHx8vPr06aNZs2bp0KFDoS6pWSoqKiRJMTExkqSCggKdOHGiwevZt29fde/e3dOv5+nnecprr72m2NhYDRgwQLm5uTp27FgoyguKkydP6o033lBVVZUyMjJa7WvpZa67KcfpDh48qJMnTyohIaHB4wkJCfryyy9DVFXwpaena+nSperTp4/279+vhx56SNdff7127dqlzp07h7q8oCsrK5OkRl/XU99rLcaMGaOJEycqNTVVu3fv1t/+7d9q7Nix2rx5s9q0aRPq8qzV1dVp9uzZuvbaazVgwABJ372eERER6tKlS4NYL7+ejZ2nJP3kJz9Rjx49lJycrB07dui+++5TUVGR3nnnnRBWa2/nzp3KyMhQdXW1OnXqpBUrVqh///4qLCxsda+l17m+UfvF2LFj6/88aNAgpaenq0ePHnrrrbc0ffr0EFaG5rrtttvq/zxw4EANGjRIvXr10oYNGzRy5MgQVtY0WVlZ2rVrl+evoTifs53nzJkz6/88cOBAJSUlaeTIkdq9e7d69ep1octssj59+qiwsFAVFRV6++23NWXKFOXn54e6LDTC9UvfsbGxatOmzRlXHJaXlysxMTFEVbW8Ll266LLLLlNxcXGoS2kRp147v72uktSzZ0/FxsZ68rXNzs7W6tWr9dFHHzW4HW1iYqJqa2t15MiRBvFefT3Pdp6NSU9PlyTPvZ4RERHq3bu30tLSlJeXp8GDB+u5555rda9la+D6Rh0REaG0tDStW7eu/rG6ujqtW7dOGRkZIaysZR09elS7d+9WUlJSqEtpEampqUpMTGzwugYCAW3durVVv66StHfvXh06dMhTr63jOMrOztaKFSu0fv16paamNvh+Wlqa2rZt2+D1LCoqUklJiadez/OdZ2MKCwslyVOvZ2Pq6upUU1PTal7LViXUV7OZeOONN5zIyEhn6dKlzueff+7MnDnT6dKli1NWVhbq0oJmzpw5zoYNG5w9e/Y4//7v/+5kZmY6sbGxzoEDB0JdWpNVVlY627dvd7Zv3+5Icp5++mln+/btzn//9387juM4jz32mNOlSxdn5cqVzo4dO5zx48c7qampzvHjx0NcuZ1znWdlZaXz61//2tm8ebOzZ88e58MPP3T++q//2rn00kud6urqUJdubNasWU50dLSzYcMGZ//+/fXHsWPH6mPuuusup3v37s769eudbdu2ORkZGU5GRkYIq7Z3vvMsLi52Hn74YWfbtm3Onj17nJUrVzo9e/Z0hg4dGuLK7cydO9fJz8939uzZ4+zYscOZO3euExYW5vzxj390HKd1vJatiScateM4zgsvvOB0797diYiIcIYMGeJs2bIl1CUF1eTJk52kpCQnIiLC+au/+itn8uTJTnFxcajLapaPPvrIkXTGMWXKFMdxvvuI1gMPPOAkJCQ4kZGRzsiRI52ioqLQFt0E5zrPY8eOOaNGjXLi4uKctm3bOj169HBmzJjhuX9kNnZ+kpwlS5bUxxw/ftz5xS9+4Vx88cVOhw4dnB/96EfO/v37Q1d0E5zvPEtKSpyhQ4c6MTExTmRkpNO7d2/n3nvvdSoqKkJbuKU777zT6dGjhxMREeHExcU5I0eOrG/SjtM6XsvWhNtcAgDgYq5/jxoAAD+jUQMA4GI0agAAXIxGDQCAi9GoAQBwMRo1AAAuRqMGAMDFaNQAALgYjRoAABejUQMA4GI0agAAXIxGDQCAi/1/9+7pn+6Nl5YAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "visualize_Q(QUBO)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "EzZFQ_MBu7Ul"
      },
      "source": [
        "各要素が正負の値が入り乱れているのが見て取れるのがわかるかと思います．それでは何か$\\vec{x}$を入力したときにこのQUBOを用いて$y$を出力するような関数を以下のように定義しておきましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "y2XtRA7wvRij"
      },
      "outputs": [],
      "source": [
        "def black_box(x):\n",
        "  assert QUBO.shape == (len(x), len(x)), \"サイズが一致しません\"\n",
        "\n",
        "  y = x @ QUBO @ x\n",
        "  return y"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uwi9Fpx-xU9k"
      },
      "source": [
        "**この関数の中身はあくまで知らないという仮定で最適化を行なっていきます．**このblack_box関数の中身を自分の好きな問題設定でやることもできるので興味がある方はぜひやってみてください．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "NBIT7QHGuOxj"
      },
      "source": [
        "## ブラックボックス最適化の流れ"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "YVAu2l1-q9KL"
      },
      "source": [
        "ブラックボックス最適化は以下のような流れで行います．\n",
        "1.   データセット$D$に最もフィットするような代理モデル$\\hat{\\; f \\;}$を作成する\n",
        "2.   $\\vec{x}^{next} = \\text{argmin}_{\\;\\;\\vec{x}} \\, \\hat{\\; f \\;}(\\vec{x})$となるようなデータをサンプリングする\n",
        "3.   $y^{next}=f(\\vec{x}^{next})$を計算する．（$f$はブラックボックス関数であるという想定）\n",
        "4.   データセット$D$に$(\\vec{x}^{next}, y^{next})$を追加する\n",
        "5.   1.-4.を繰り返す"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "1dUxQWbmZIyq"
      },
      "source": [
        "### 1. 代理モデル$\\hat{\\; f \\;}$をデータセット$D$を用いて学習を行う"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Ee5pY4hgaeFc"
      },
      "source": [
        "実際に定式化ができない問題がある場合でもとにかく何かしら代わりの式を用意しないことには何もできません．このような式を代理モデルといいます．そして今回，代理モデルとして以下のような2次の多項式を用います．\n",
        "\\begin{equation}\n",
        "  \\hat{\\; f \\;}(\\vec{x}) = a_0 + \\sum_{i=1}^N a_i x_i + \\sum_{i=1}^N \\sum_{j=i+1}^N a_{ij} x_i x_j\n",
        "\\end{equation}\n",
        "この式はまさにQUBOの形です．つまり，量子アニーリングなどでこの代理モデルの最適化を行うことができます．  \n",
        "そして，この式をより簡単な式に変形させましょう．そのために以下のような変換を行います．\n",
        "\\begin{equation}\n",
        "  \\vec{z} = (1, x_1, x_2, ..., x_N, x_1x_2, x_1x_3, ..., x_{N-1}x_N)^{\\rm T} \\\\\n",
        "  \\vec{a} = (a_0, a_1, a_2, ..., a_N, a_{12}\\ ,a_{13}\\ , ..., a_{N-1N})^{\\rm T}\n",
        "\\end{equation}\n",
        "すると，代理モデル$\\hat{\\; f \\;}(\\vec{x})$は以下のように変形することができます．\n",
        "\\begin{eqnarray*}\n",
        "  \\hat{\\; f \\;}(\\vec{x}) &=& a_0 + \\sum_{i=1}^N a_i x_i + \\sum_{i=1}^N \\sum_{j=i+1}^N a_{ij} x_i x_j \\\\\n",
        "  &=& a + a_1x_1 + a_2x_2 + \\ldots + a_Nx_N + a_{12}\\ x_1x_2 + a_{13}\\ x_1x_3 + \\ldots + a_{N-1N}x_{N-1}x_N \\\\\n",
        "  &=& \\vec{a}^{\\rm T} \\vec{z}\n",
        "\\end{eqnarray*}\n",
        "これで代理モデルを線形モデルとして記述することができました．線形モデルで記述できるメリットは線形回帰を用いることができる点です．線形回帰とはいくつかの入力$\\vec{z}$に対して一次関数の範囲で最も合うような$\\vec{a}$を見つける手法です．"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dRnIvHe76N_8"
      },
      "source": [
        "ここから，実際にコードで記述してみましょう．まずは，$\\vec{z}$の変数の数は\n",
        "\\begin{eqnarray}\n",
        "  P &=& 1 + N + {}_N C_2 \\\\\n",
        "  &=& 1 + N + \\frac{N(N-1)}{2}\n",
        "\\end{eqnarray}\n",
        "となります．まずはこれを計算するための関数を定義しましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "xLbZa9TfA3Nw"
      },
      "outputs": [],
      "source": [
        "def calc_P_from_N(N):\n",
        "  P = 1 + N + N * (N-1) / 2\n",
        "  return int(P)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "Zrp6A4ftBBR0"
      },
      "source": [
        "次に$\\vec{x}$から$\\vec{z}$に変換を行う関数と$\\vec{a}$からQUBOに変換する関数は以下のようになります.\n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "L0Dor5mj5cnU"
      },
      "outputs": [],
      "source": [
        "def calc_z_from_x(x):\n",
        "  N = len(x) # 変数の数\n",
        "  P = calc_P_from_N(N) # Xの変数の数\n",
        "  z = np.zeros(P)\n",
        "\n",
        "  z[0] = 1\n",
        "  z[1:N+1] = x\n",
        "\n",
        "  n = 0\n",
        "  for i in range(N-1):\n",
        "    for j in range(i+1, N):\n",
        "      z[N+1+n] = x[i] * x[j]\n",
        "      n = n + 1\n",
        "\n",
        "  return z\n",
        "\n",
        "def calc_Q_from_a(a):\n",
        "  # 対角成分にa_1~a_Nを代入する\n",
        "  Q = np.diag(a[1:N+1])\n",
        "\n",
        "  n = 0\n",
        "  for i in range(N-1):\n",
        "    for j in range(i+1, N):\n",
        "      Q[i,j] = a[N+1+n]\n",
        "      n = n + 1\n",
        "\n",
        "  return Q"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "BDXz0kgMn7aL"
      },
      "source": [
        "$\\vec{a}$を線形回帰で求める場合，以下のように表されます．  \n",
        "\\begin{equation}\n",
        "  \\vec{a} = {\\rm arg}\\min_{\\vec{a}} \\left(\n",
        "    \\sum_{i=1}^n (y^{(i)} - \\vec{a}^{\\rm T} \\vec{z}^{(i)})^2 + \\lambda \\vec{a}^{\\rm T}\\vec{a}\n",
        "    \\right)\n",
        "\\end{equation}\n",
        "ここで，$n$はデータの数を表しています．第1項が真のモデルと代理モデルの差を最小にするための項です．これは最小二乗法に当たります．そして，第2項が正則化項と呼ばれるものです．正則化項は過学習を防ぐために用いられます．  \n",
        "それでは与えられたデータセット$D$から$\\vec{a}$を推測する関数を定義しましょう．  \n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "pfVI3tWkh2Gm"
      },
      "outputs": [],
      "source": [
        "def predict_a(z_data, y_data, alpha):\n",
        "  \"\"\"\"\n",
        "  Args:\n",
        "    z_data: D個のzのデータ (D, P)\n",
        "    y_data: D個のyのデータ (D, )\n",
        "  Return:\n",
        "    a: 回帰係数 (P, )\n",
        "  \"\"\"\n",
        "  # 線形回帰\n",
        "  ridge = Ridge(alpha=alpha)\n",
        "  ridge.fit(z_data, y_data)\n",
        "  a = ridge.coef_\n",
        "\n",
        "\n",
        "  return a"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "dvKxlp_4xZMz"
      },
      "source": [
        "### 2. $\\vec{x}^{next} = \\text{argmin}_{\\;\\;\\vec{x}} \\, \\hat{\\; f \\;}(\\vec{x})$となるようなデータをサンプリングする"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "7u3o-Q9jxeZM"
      },
      "source": [
        "1.で今持っているデータセットをもとに線形回帰を行うことで代理モデルを作りました．次にすべきことは作成した代理モデルの最適解を求めることです．今回代理モデルとしておいた形はQUBOと同じだと先ほど述べました．ですので，QUBOを最適化する手法としてSAやQAを用いましょう．  \n"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "JREAk594jDzw"
      },
      "outputs": [],
      "source": [
        "def sample_from_qubo(sampler, qubo):\n",
        "  \"\"\"\n",
        "  Args:\n",
        "    sampler: サンプラー（SAやQAなど）\n",
        "    qubo   : QUBO (N, N)\n",
        "  Return:\n",
        "    x      : 最適解 (N, )\n",
        "  \"\"\"\n",
        "  sampleset = sampler.sample_qubo(qubo, num_reads=100)\n",
        "  # 1番エネルギーが低いサンプル\n",
        "  sample = sampleset.first.sample\n",
        "\n",
        "  # 出力は辞書形式なのでそれをリストに変換\n",
        "  x = []\n",
        "  for i in range(len(sample)):\n",
        "    x.append(sample[i])\n",
        "\n",
        "  return x"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### ブラックボックス最適化の実装"
      ],
      "metadata": {
        "id": "EeUIKNDPU9WD"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "J5ce0-Xl5Gt2"
      },
      "source": [
        "ここまできたら，あとは3. で最初に定義したblack_box関数に2.で求めた最適解$\\vec{x}^{next}$を代入して出てきた$y^{next}$を求めます．最後に4. で$(\\vec{x}^{next}, y^{next})$をデータセットとして元々持っているデータセットに加えていくだけです．  \n",
        "なのであとは1.~4.を繰り返しましょう．  \n",
        "もう一度，ブラックボックス最適化の流れを示します．  \n",
        "1.   データセット$D$に最もフィットするような代理モデル$\\hat{\\; f \\;}$を作成する\n",
        "2.   $\\vec{x}^{next} = \\text{argmin}_{\\;\\;\\vec{x}} \\, \\hat{\\; f \\;}(\\vec{x})$となるようなデータをサンプリングする\n",
        "3.   $y^{next}=f(\\vec{x}^{next})$を計算する．（$f$はブラックボックス関数であるという想定）\n",
        "4.   データセット$D$に$(\\vec{x}^{next}, y^{next})$を追加する\n",
        "5.   1.-4.を繰り返す"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "E3I-gt8D8uKq"
      },
      "source": [
        "まずは最初は適当でいいのでデータセットを用意する必要があります．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "TCoCutXXjP-u"
      },
      "outputs": [],
      "source": [
        "init_data_num = 5 # 初期データセット数\n",
        "T_all = 256 # 反復回数\n",
        "sampler = SimulatedAnnealingSampler()\n",
        "\n",
        "# データセットを作成\n",
        "init_x_data = np.random.randint(0, 2, (init_data_num, N))\n",
        "init_y_data = np.zeros(init_data_num)\n",
        "for i in range(init_data_num):\n",
        "  init_y_data[i] = black_box(init_x_data[i])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "gTiO586B25Yl"
      },
      "outputs": [],
      "source": [
        "# x_dataとy_dataをリストに変換\n",
        "x_data = init_x_data.tolist()\n",
        "y_data = init_y_data.tolist()\n",
        "\n",
        "# zを計算しておく\n",
        "z_data = []\n",
        "for i in range(len(x_data)):\n",
        "  z = calc_z_from_x(x_data[i])\n",
        "  z_data.append(z)\n",
        "\n",
        "# yの中で最もエネルギーが低い値を保存\n",
        "y_data.sort(reverse=True)\n",
        "y_min = y_data[-1]\n",
        "y_best_data = []\n",
        "y_best_data.append(y_min)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TjCqY2iw849m"
      },
      "source": [
        "反復アルゴリズムは以下のようになります．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "P8Xyg0Vy7Lq5",
        "outputId": "e84e5670-bf4c-411d-b471-e1291a0854e7"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "100%|██████████| 256/256 [00:51<00:00,  4.96it/s]\n"
          ]
        }
      ],
      "source": [
        "for t in tqdm(range(T_all)):\n",
        "  # 1. 代理モデルを作成する\n",
        "  # データセットから最もフィットする係数ベクトルaを求める\n",
        "  a = predict_a(z_data, y_data, alpha=0.01)\n",
        "  # 係数ベクトルaからQUBOに変換する\n",
        "  Q = calc_Q_from_a(a)\n",
        "\n",
        "  # 2. QUBOを最適解を求める\n",
        "  x_next = sample_from_qubo(sampler, Q)\n",
        "\n",
        "  # 3. y_next=f(x_next)を求める\n",
        "  y_next = black_box(x_next)\n",
        "\n",
        "  # 4. データセットに追加する\n",
        "  x_data.append(x_next)\n",
        "  y_data.append(y_next)\n",
        "\n",
        "  z_next = calc_z_from_x(x_next)\n",
        "  z_data.append(z_next)\n",
        "\n",
        "  # 現時点のデータセットで最もエネルギーが低いyをy_best_dataに追加\n",
        "  if y_next < y_best_data[-1]: # 一番最後尾の値と比較\n",
        "    y_best_data.append(y_next)\n",
        "  else:\n",
        "    y_best_data.append(y_best_data[-1])"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "TglAsiBvz60v"
      },
      "source": [
        "また，ブラックボックス関数であるQUBOの近似的な最適解も求めておきましょう．  \n",
        "この解からどれくらい離れているかを見ることによって，ある時点における性能を測ることができます．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "OP_Fbl393Bu_",
        "outputId": "8d0d6d91-7a93-471c-c44c-83ec8ef1c0a5"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "-32.69375799050562\n"
          ]
        }
      ],
      "source": [
        "x = sample_from_qubo(sampler, QUBO)\n",
        "optim = black_box(x)\n",
        "print(optim)"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 結果の可視化"
      ],
      "metadata": {
        "id": "EFivLXShVNzW"
      }
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "6epY78T1gHN4"
      },
      "source": [
        "それでは実際に最適解からどれくらい離れているかを見るためにグラフをプロットするための関数を定義しておきましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "vVYgRSCOQ6Uk"
      },
      "outputs": [],
      "source": [
        "def plot_progress(y_best_data, y_data, optim, file_name):\n",
        "  \"\"\"\n",
        "  Args:\n",
        "    y_best_data: 現時点における最適解のリスト\n",
        "    y_data     : 現時点におけるyのリスト\n",
        "    optim      : 最適解\n",
        "  \"\"\"\n",
        "  plt.cla()\n",
        "  plt.plot(range(len(y_best_data)), y_best_data-optim, color=\"b\", label=\"best\")\n",
        "  plt.scatter(range(len(y_data)), y_data-optim, c=\"r\", marker=\"x\", label=\"sample\")\n",
        "  # plt.ylim(top=50, bottom=min(y_best_data)-optim-1)\n",
        "  plt.ylim(top=50, bottom=-1)\n",
        "  plt.grid()\n",
        "  plt.axhline(y=0, linestyle=\"--\", c=\"g\")\n",
        "  plt.xlabel(\"step\")\n",
        "  plt.ylabel(\"optim - energy\")\n",
        "  plt.legend(bbox_to_anchor=(1, 1), loc='upper right')\n",
        "  plt.savefig(file_name)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "uSM3iG5K1C3N"
      },
      "source": [
        "また，折角ならアニメーションでブラックボックス最適化の様子を見てみたいですよね．なのでそのための準備としてフォルダを作っておきましょう．注意すべきなことは，このフォルダは仮想的なものなのでランタイムが切断されたりすると消えてしまいます．消えてしまった場合でも同じように実行するとまた新しく作られるので安心してください．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "7MyvOfoqbSwI",
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "outputId": "9469c49c-fd9c-4472-a1af-2bd98be8bc5a"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "mkdir: cannot create directory ‘fig_data’: File exists\n"
          ]
        }
      ],
      "source": [
        "!mkdir fig_data"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "U-wQHj4o1pG0"
      },
      "source": [
        "それでは先ほど作ったfig_dataに1ステップずつグラフをプロットしたものを格納していきましょう．"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 472
        },
        "id": "q-wmXwzPdcvG",
        "outputId": "4f6a6c81-c611-4f03-d57c-ae6c11c7e259"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "100%|██████████| 256/256 [01:05<00:00,  3.92it/s]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAG2CAYAAABlBWwKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCbUlEQVR4nO3de1yUdd7/8fcAw1nAEyAKqOUpUysrIzuY4qlsLb23Mje1LG8LzUO25W5Zupu03XdptVa71drh3g6rW7aVa5IFZh7ynFn58xiYHDLjJAIDXL8/ZhmdAGVwLmYuej0fDx461/Wd73zmwzC8ua5rrstmGIYhAAAACwrwdQEAAABNRZABAACWRZABAACWRZABAACWRZABAACWRZABAACWRZABAACWRZABAACWRZABAACWRZABAACW5dMg8+ijj8pms7l99ezZ07W+vLxcaWlpatu2rSIjIzV27Fjl5+f7sGIAAOBPfL5Fpnfv3srNzXV9rVu3zrVu1qxZev/997Vs2TJlZWXpyJEjGjNmjA+rBQAA/iTI5wUEBSk+Pr7O8qKiIr388st64403NHjwYEnS0qVL1atXL23cuFGXXXZZc5cKAAD8jM+DzN69e5WQkKDQ0FClpKQoPT1dSUlJ2rp1qxwOh1JTU11je/bsqaSkJG3YsKHBIFNRUaGKigrX7ZqaGh07dkxt27aVzWYz/fkAAICzZxiGSkpKlJCQoICAhncg+TTIDBgwQK+88op69Oih3NxczZ8/X1deeaW++uor5eXlKTg4WDExMW73iYuLU15eXoNzpqena/78+SZXDgAAmkNOTo46derU4HqfBpmRI0e6/t+3b18NGDBAycnJ+sc//qGwsLAmzTl37lzNnj3bdbuoqEhJSUk6ePCgWrVqddY113I4HPr00091zTXXyG63e21enESPzUV/zUePzUV/zefLHpeUlKhLly5n/N3t811Lp4qJiVH37t21b98+DR06VJWVlSosLHTbKpOfn1/vMTW1QkJCFBISUmd5mzZtFBUV5bVaHQ6HwsPD1bZtW36ATEKPzUV/zUePzUV/zefLHtc+3pkOC/H5p5ZOVVpaqv3796tDhw7q37+/7Ha71qxZ41q/Z88eZWdnKyUlxYdVAgAAf+HTLTJz5szR9ddfr+TkZB05ckSPPPKIAgMDNW7cOEVHR2vy5MmaPXu2a2vK9OnTlZKSwieWAACAJB8HmcOHD2vcuHH68ccf1b59e11xxRXauHGj2rdvL0latGiRAgICNHbsWFVUVGj48OF67rnnfFkyAADwIz4NMm+99dZp14eGhmrJkiVasmRJM1UEAPglqK6ulsPh8HUZfs/hcCgoKEjl5eWqrq726tx2u12BgYFnPY9fHewLAICZDMNQXl6eCgsLfV2KJRiGofj4eOXk5JhyLraYmBjFx8ef1dwEGQDAL0ZtiImNjVV4eDgnSj2DmpoalZaWKjIy8rQnpfOUYRgqKytTQUGBJKlDhw5NnosgAwD4RaiurnaFmLZt2/q6HEuoqalRZWWlQkNDvRpkJLnOF1dQUKDY2Ngm72byq49fAwBgltpjYsLDw31cCWrVfi/O5nglggwA4BeF3Un+wxvfC4IMAACwLIIMAAB+btCgQZo5c6avy/BLBBkAAH7BMjMzZbPZLPuRdIIMAACwLIIMAAAWUFVVpWnTpik6Olrt2rXTww8/LMMwJEkVFRWaM2eOOnbsqIiICA0YMECZmZmu+3733Xe6/vrr1bp1a0VERKh3795auXKlDh06pGuuuUaS1Lp1a9lsNk2aNMkHz67pOI8MAOAXyzCksjLfPHZ4uOTJh3ZeffVVTZ48WV988YW2bNmiKVOmKCkpSXfddZemTZumr7/+Wm+99ZYSEhL07rvvasSIEdq1a5e6deumtLQ0VVZWau3atYqIiNDXX3+tyMhIJSYm6p///KfGjh2rPXv2KCoqynV+F6sgyAAAfrHKyqTISN88dmmpFBHR+PGJiYlatGiRbDabevTooV27dmnRokUaPny4li5dquzsbCUkJEiS5syZo1WrVmnp0qVauHChsrOzNXbsWPXp00eS1LVrV9e8bdq0kSTFxsYqJibGa8+vuRBkAACwgMsuu8ztvCspKSl68skntWvXLlVXV6t79+5u4ysqKlxnML733nt19913a/Xq1UpNTdXYsWPVt2/fZq3fLAQZAMAvVni4c8uIrx7bG0pLSxUYGKitW7fWOc1/5H82N915550aPny4PvzwQ61evVrp6el68sknNX36dO8U4UMEGQDAL5bN5tnuHV/atGmT2+2NGzeqW7duuvDCC1VdXa2CggJdeeWVDd4/MTFRU6dO1dSpUzV37ly9+OKLmj59uoKDgyU5r0VlRXxqCQAAC8jOztbs2bO1Z88evfnmm3r22Wc1Y8YMde/eXePHj9eECRP0zjvv6ODBg/riiy+Unp6uDz/8UJI0c+ZMffTRRzp48KC2bdumTz/9VL169ZIkJScny2az6YMPPtAPP/ygUl9tomoigow3FBVJhw/Xv+7wYed6AADOwoQJE3TixAldeumlSktL04wZMzRlyhRJ0tKlSzVhwgTdd9996tGjh2644QZt3rxZSUlJkpxbW9LS0tSrVy+NGDFC3bt313PPPSdJ6tixo+bPn68HH3xQcXFxmjZtms+eY1Owa+lsFRdLo0ZJBQVSZqaUmHhyXU6ONGiQFBsrrVolRUf7qkoAgIWdek6Y559/vs56u92u+fPna/78+fXe/9lnnz3t/A8//LAefvjhs6rRV9gic7ZKS50h5sABZ2jJyXEurw0xBw4415eU+LJKAABaJILM2UpIcG6J6dr1ZJhZv/5kiOna1bm+Uyff1gkAQAvEriVvSEx0hpXa8DJwoHN5bYg5dXcTAADwGrbIeEtiovT66+7LXn+dEAMAgIkIMt6SkyPddpv7sttuO3nMDAAA8DqCjDecemBv167S55+7HzNDmAEAwBQEmbN15EjdA3svv7zuAcANnWcGAAA0GQf7nq3ISOd5YiT3A3tPPQA4NlZq1cpHBQIA0HIRZM5WVJTzZHclJXU/Yp2YKGVlOUMMJ8MDAMDrCDLeEB3dcFDh/DEAgBZq0qRJKiws1IoVK3xWA8fIAAAAyyLIAADQWFwk2O8QZAAAaIyiImnECOnqq+ueViMnx7l8xAjTwszy5cvVp08fhYWFqW3btkpNTdXx48e1efNmDR06VO3atVN0dLSuvvpqbdu2ze2+NptNf/nLXzRq1CiFh4erV69e2rBhg/bt26dBgwYpIiJCl19+ufbv3++6z6OPPqqLLrpIS5cuVXJyssLDw3XTTTep6DTPr6amRunp6erSpYvCwsLUr18/LV++3JR+1CLIeBNJHQBarpISn10kODc3V+PGjdMdd9yhb775RpmZmRozZowMw1BJSYkmTpyodevWaePGjerWrZuuvfZalfysjj/84Q+aMGGCduzYoZ49e+rWW2/Vf//3f2vu3LnasmWLDMPQtGnT3O6zb98+rVixQu+9955WrVql7du365577mmwzvT0dL322mt64YUXtHv3bs2aNUu/+c1vlJWV5fWe1OJgX2+pTeoFBXWvr1T7Io+NdX7CiU8wAYD1dOrkfl29QYOcl6K57TbTLxKcm5urqqoqjRkzRsnJyZKkPn36SJIGDx7sNvavf/2rYmJilJWVpVGjRrmW33777brpppskSQ888IBSUlL08MMPa/jw4ZKkGTNm6Pbbb3ebq7y8XM8//7x69uypgIAAPfvss7ruuuv05JNPKj4+3m1sRUWFFi5cqI8//lgpKSmSpK5du2rdunX6y1/+oquvvtqLHTmJLTLe4sOkDgBoJrXnCKs94enAge4hxqTr6/Xr109DhgxRnz599Otf/1ovvviifvrpJ0lSfn6+7rrrLnXr1k3R0dGKiopSaWmpsrOz3ebo27ev6/9xcXGSToah2mXl5eUqLi52LUtKSlJCQoLrdkpKimpqarRnz546Ne7bt09lZWUaOnSoIiMjXV+vvfaa2y4rb2OLjLf4MKkDAJpR7UWCBw48uczkiwQHBgYqIyND69ev1+rVq/Xss8/q97//vTZt2qS7775bP/74o55++mklJycrJCREKSkpqqysdJvDbre7/m+z2RpcVlNT06QaS0tLJUkffvihOnbs6LYuJCSkSXM2BkHGm049m29tUpdMT+oAgGbU0EWCTX6ft9lsGjhwoAYOHKh58+YpOTlZ7777rj7//HM999xzuvbaa/9TXo6OHj3qlcfMzs5Wbm6uoqKiJEkbN25UQECAevToUWfseeedp5CQEGVnZ5u2G6k+BBlv80FSBwA0k59fJPjULe+DBpkWZjZt2qQ1a9Zo2LBhio2N1aZNm/TDDz+oV69e6tatm15//XVdfPHFKi4u1v3336+wsDCvPG5oaKjuuecePfXUUyotLdW9996rm266qc7xMZLUqlUrzZkzR7NmzVJNTY2uuOIKFRUV6fPPP1dUVJQmTpzolZp+jiDjbT5K6gAAkx0+XPciwT/fEj9okPPSNF4+jCAqKkpr167V4sWLVVxcrOTkZD355JMaOXKk4uPjNWXKFF100UVKTEzUwoULNWfOHK887rnnnqtRo0Zp1KhROnbsmEaNGqXnnnuuwfF/+MMf1L59e6Wnp+vAgQOKiYnRRRddpN/97ndeqac+BBlv8lFSBwA0g1atfHaR4F69emnVqlX1rrvwwgu1efNmt2X/9V//5XbbMAy32507d66zbNCgQXWWSdLkyZM1a9YsBQTU/XzQK6+84nbbZrNpxowZmjFjRoPPxdsIMt7iw6QOAGgG0dFcJNgPEWS8xYdJHQDQTLhIsN/hPDLeUpvUs7Lq7j6qTeqcDA8AYBGPPvponUsd+CO2yHgTSR0AgGbFFhkAwC9KfQe0wje88b0gyAAAfhFqz2JbVlbm40pQq/Z7ceoZhj3FriUAwC9CYGCgYmJiVFBQIEkKDw93nZYf9aupqVFlZaXKy8vr/fh1UxmGobKyMhUUFCgmJkaBgYFNnosgAwD4xag9I21tmMHpGYahEydOKCwszJTQFxMTU+9Zgj1BkAEA/GLYbDZ16NBBsbGxcjgcvi7H7zkcDq1du1ZXXXXVWe3+qY/dbj+rLTG1CDIAgF+cwMBAr/wSbekCAwNVVVWl0NBQrwcZb+FgXwAAYFkEGQAAYFkEGQAAYFkEGQAAYFkEGQAAYFkEGQAAYFkEGQAAYFkEGQAAYFkEGQAAYFkEGQAAYFkEGQAAYFkEGQAAYFl+E2Qef/xx2Ww2zZw507WsvLxcaWlpatu2rSIjIzV27Fjl5+f7rkgAAOBX/CLIbN68WX/5y1/Ut29ft+WzZs3S+++/r2XLlikrK0tHjhzRmDFjfFQlAADwNz4PMqWlpRo/frxefPFFtW7d2rW8qKhIL7/8sp566ikNHjxY/fv319KlS7V+/Xpt3LjRhxUDAAB/4fMgk5aWpuuuu06pqaluy7du3SqHw+G2vGfPnkpKStKGDRuau0wAAOCHgnz54G+99Za2bdumzZs311mXl5en4OBgxcTEuC2Pi4tTXl5eg3NWVFSooqLCdbu4uFiS5HA45HA4vFP4f+Y79V94Hz02F/01Hz02F/01ny973NjH9FmQycnJ0YwZM5SRkaHQ0FCvzZuenq758+fXWb569WqFh4d77XFqZWRkeH1OuKPH5qK/5qPH5qK/5vNFj8vKyho1zmYYhmFyLfVasWKFbrzxRgUGBrqWVVdXy2azKSAgQB999JFSU1P1008/uW2VSU5O1syZMzVr1qx6561vi0xiYqKOHj2qqKgor9XvcDiUkZGhoUOHym63e21enESPzUV/zUePzUV/zefLHhcXF6tdu3YqKio67e9vn22RGTJkiHbt2uW27Pbbb1fPnj31wAMPKDExUXa7XWvWrNHYsWMlSXv27FF2drZSUlIanDckJEQhISF1ltvtdlO+CWbNi5Posbnor/nosbnor/l80ePGPp7PgkyrVq10/vnnuy2LiIhQ27ZtXcsnT56s2bNnq02bNoqKitL06dOVkpKiyy67zBclAwAAP+PTg33PZNGiRQoICNDYsWNVUVGh4cOH67nnnvN1WY1TVCSVlEidOtVdd/iw1KqVFB3d/HUBANCC+FWQyczMdLsdGhqqJUuWaMmSJb4pqDGOHJGSk92XFRVJgwdLx45Ja9dKiYkn1+XkSIMGSbGx0qpVhBkAAM6Cz88jY1n/+Vi3Ro50hpNT7dkjffmldOiQdNVVJ9fXhpgDB6SCAucWGwAA0GQEmaYqLXX+e+iQM5zUhpXdu6Vf/1qqqpKCgk6uX79euuIKZ4jp2lXKzKx/txMAAGg0gkxTJSQ4/+3c2RlOBg2SVq+W+veXsrOlpCRp3TpnaDlwQBo40Lk8JET617/cdzcBAIAmIcg0RVGR89gYSVq58mRYGT5cOuUcNkpIkF5/3f2+bdtyXAwAAF5CkPFUUZE0YoTz2BhJ6tixblix251bX664wrmb6VRBQZJvzkEIAECLQ5DxVEmJ80DdQ4ect7dskW67zX1MTY0UF+cMM7VbbhISnLubsrPdj6kBAABNRpDxVKdOzgN1O3d23h4+3LlbKeg/n2QPCpKqq6WjR93vt2yZ+zEzgwY5zycDAACajCDTFImJ0t/+5vx/VdXJf7t2dYaVTp2cYeZUtVttMjOd42JjnSfFAwAATUaQaapu3eoue/115y6kgP+0NThY+ugj960wkpSVxcnwAADwAoJMU9V3Mrtx45wH+NZ+/HrbNmnYsJNbYU4NM4QYAADOGkGmKXJypGuvdf6/c2fp88+dQSU7W8rPP3kOmd69nWMSE9mlBACACfzqWkuWcPiwc6tKbq7z9sqVzjCTmXny8gOSZLO53y8x0blLiYtFAgDgNQQZT7Vq5dyqEhLivN2xo/Pf2q0utReErG+rC5ckAADAqwgynoqOdh6oW1go7djhvo6tLgAANCuOkWmK6OiT11r6uU6dCDEAADQTggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsggwAALAsgownioqkw4frX3f4sHM9AABoNgSZxioqkkaMkK6+WsrJcV+Xk+NcPmIEYQYAgGZEkGmskhKpoEA6cEAaNEj6/nvn8u+/d94+cMC5vqTEl1UCAPCLQpBprE6dpMxMqWtXZ2i59lrn8muvdd7u2tW5vlMnX1YJAMAvSpCvC7CUxERnWBk0SI/snaBX7h6iiLx/SUFBUk1HabC9wbumpEivvirZbM1XLgAALR1BxlOJidLrr+voVV8rNzdS0rlSlaRDp7/bvn3SwoVssAEAwJsIMp7KyZFuu02zg0J1zmPdlTLvEQW1ayf9+c9SXFy9dxkxwnnoTGlpM9cKAEALR5DxRE6O68Dec3r1Us+ePynlnHzZv8mU5mxz7nZKTKxzt6goZ5A5fry5CwYAoGXjYN/GOnz45KeTunaVVq50Ll+58uQBwIMG1XuemYgI578EGQAAvIsg01itWkmxsSc/ndSxo3N5x44nP80UG+sc9zMEGQAAzMGupcaKjpZWrXLuI+rUSXI4Tq5LTJSyspwhJjq6zl0JMgAAmIMg44no6HqDiqTTfhyJIAMAgDnYtdQMIiOd/xJkAADwLoJMM6jdIsPHrwEA8C6CTDNg1xIAAOYgyDQDggwAAOYgyDQDggwAAObwOMhMnDhRa9euNaOWFosgAwCAOTwOMkVFRUpNTVW3bt20cOFCff/992bU1aIQZAAAMIfHQWbFihX6/vvvdffdd+vtt99W586dNXLkSC1fvlyOU08SBxeCDAAA5mjSMTLt27fX7NmztXPnTm3atEnnnnuubrvtNiUkJGjWrFnau3evt+u0NIIMAADmOKuDfXNzc5WRkaGMjAwFBgbq2muv1a5du3Teeedp0aJF3qrR8ggyAACYw+Mg43A49M9//lOjRo1ScnKyli1bppkzZ+rIkSN69dVX9fHHH+sf//iHFixYYEa9lkSQAQDAHB4HmQ4dOuiuu+5ScnKyvvjiC23ZskVTp05VVFSUa8w111yjmJiYM871/PPPq2/fvoqKilJUVJRSUlL073//27W+vLxcaWlpatu2rSIjIzV27Fjl5+d7WrLPEWQAADCHx0Fm0aJFOnLkiJYsWaILLrig3jExMTE6ePDgGefq1KmTHn/8cW3dulVbtmzR4MGDNXr0aO3evVuSNGvWLL3//vtatmyZsrKydOTIEY0ZM8bTkn2OIAMAgDk8vvr1bbfd5rUHv/76691uP/bYY3r++ee1ceNGderUSS+//LLeeOMNDR48WJK0dOlS9erVSxs3btRll13mtTrMVnvRSK61BACAd3kcZBraImKz2RQaGqpzzz1Xt956q3r06OHRvNXV1Vq2bJmOHz+ulJQUbd26VQ6HQ6mpqa4xPXv2VFJSkjZs2NBgkKmoqFBFRYXrdnFxsSTnsT3e/Hh47VyNmTM4WJLsqqqSjh93/Oc2zsSTHsNz9Nd89Nhc9Nd8vuxxYx/T4yATFRWlFStWKCYmRv3795ckbdu2TYWFhRo2bJjefvtt/elPf9KaNWs0cODAM863a9cupaSkqLy8XJGRkXr33Xd13nnnaceOHQoODq5zrE1cXJzy8vIanC89PV3z58+vs3z16tUKDw/37Mk2QkZGxhnHOBw2Sb+SJK1YkaHISH7oPNGYHqPp6K/56LG56K/5fNHjsrKyRo3zOMjEx8fr1ltv1Z///GcFBDgPsampqdGMGTPUqlUrvfXWW5o6daoeeOABrVu37ozz9ejRQzt27FBRUZGWL1+uiRMnKisry9OyXObOnavZs2e7bhcXFysxMVHDhg1zOyD5bDkcDmVkZGjo0KGy2+1nHB8UZKiqyqbLLx+qTp28VkaL5mmP4Rn6az56bC76az5f9rh2j8qZeBxkXn75ZX3++eeuECNJAQEBmj59ui6//HItXLhQ06ZN05VXXtmo+YKDg3XuuedKkvr376/Nmzfr6aef1s0336zKykoVFha6bZXJz89XfHx8g/OFhIQoJCSkznK73W7KN6Gx80ZESEVFUmWlXfy8ecas7x2c6K/56LG56K/5fNHjxj6ex59aqqqq0rfffltn+bfffqvq6mpJUmhoqGw2m6dTS3Ju3amoqFD//v1lt9u1Zs0a17o9e/YoOztbKSkpTZrbl/jkEgAA3tekTy1NnjxZv/vd73TJJZdIkjZv3qyFCxdqwoQJkqSsrCz17t37jHPNnTtXI0eOVFJSkkpKSvTGG28oMzNTH330kaKjozV58mTNnj1bbdq0UVRUlKZPn66UlBRLfWKpFkEGAADv8zjILFq0SHFxcXriiSdcJ6eLi4vTrFmz9MADD0iShg0bphEjRpxxroKCAk2YMEG5ubmKjo5W37599dFHH2no0KGuxwoICNDYsWNVUVGh4cOH67nnnvO0ZL9AkAEAwPs8CjJVVVV64403dOedd+r3v/+960Ccnx9Em5SU1Kj5Xn755dOuDw0N1ZIlS7RkyRJPyvRLBBkAALzPo2NkgoKCNHXqVJWXl0uS69ICODOCDAAA3ufxwb6XXnqptm/fbkYtLRpBBgAA7/P4GJl77rlH9913nw4fPqz+/fsrovY39H/07dvXa8W1JAQZAAC8z+Mgc8stt0iS7r33Xtcym80mwzBks9lcH8GGu9rrLRFkAADwHo+DTGOuao26arfIcOFIAAC8x+Mgk5ycbEYdLR67lgAA8D6PD/aVpNdff10DBw5UQkKCvvvuO0nS4sWL9d5773m1uJaEIAMAgPd5vEXm+eef17x58zRz5kw99thjrmNiYmJitHjxYo0ePdrrRVpaUZFUUqKIgLaSwnT8aJlU5JBKSqRTL4i1b5/0n2tO1ZGX51y3f78UF1f/mNPdv6nrzHzcxsy9d6/z32+/lQIDrVO3VeauPZ6ttr8t8TXk67l/3mOr1N3cczf1cWv7u3Kl1K2bd+c2s24rzV3fa7j2/lFRUqtWzt9ntf8WF0sdO0rR0fXPbwKPg8yzzz6rF198UTfccIMef/xx1/KLL75Yc+bM8WpxlldUJI0YIeXlKaL6NkkLdHzNJmnwHKmgwPllGM6vqiopKEgKCHDettmc/0qSwyGFhEgVFfWPOd39m7rOzMdt7Nx2u/Tmm9KVV0plZdap2ypzh4ZKr70mXXGFs78t8TXk67kjIqRXXz35GrZK3VZ5H6p9DY8b5xznL8+pJc3989dw7f3tdmfQad9e+uknqXVr6YcfnF99+0oZGc0WZpp0sO+FF15YZ3lISIiOs9/EXUmJM6wcOqSIgD2SpOMVgdLOnSdT7qmqqhqeq6LizGPMWGfm455p7torn1ZWOn+4vDn32axrKXPX/nXlcJxcZvXn5G9zV1Y6///z17C/1+2LuZvyuKduqfXH59QS5m7oNexwSIcPS7m5zt9nOTknf6/l5zt//zVTkPH4GJkuXbpox44ddZavWrVKvXr18kZNLUenTlJmptS1qyJqSiRJxxVRf4gBAMBqqqudgbL291pSkrRunfP3XzPxeIvM7NmzlZaWpvLychmGoS+++EJvvvmm0tPT9dJLL5lRo7UlJkqZmYoY8FspV9qm/gpTma+rsoYTUs1NwQqo/NHXlbRM9Nd89Nhc9Nd8jelxtfRnTdPkpI+dISYxsfnqUxOCzJ133qmwsDA99NBDKisr06233qqEhAQ9/fTTrpPl4WcSE9X7pVmKvq5QRYpRucJ8XZF1VEqiX+ahv+ajx+aiv+ZrRI+rFeg8prGZQ4zUhCAjSePHj9f48eNVVlam0tJSxcbGeruuliUnR+2mj9P3ytNRtfN1NZbhCAnVp88+q2umT5e9otzX5bQ49Nd89Nhc9Nd8je1xGx2Txn1kjS0ypwoPD1d4eLi3ammZcnKkQYOkAwcUERioiOpsX1dkGY6AMMXGnlByQI7sOuHrcloc+ms+emwu+mu+RvW49hiZ7FLnpyCbOcx4fLBvfn6+brvtNiUkJCgoKEiBgYFuXzjF4cOuEON2MBR9AgC0BLW/22p/r2VnO8PM4cPNVoLHW2QmTZqk7OxsPfzww+rQoYNsNpsZdbUMrVpJsbFSTY3zY2jffuv8zH1sLOeRaczctYKDnXNZpW6rzF378Xa7/eRHsK3+nPxt7uBg99ewVeq2yvtQ7WtY8q/n1JLm/vlruPb+pzuPTFyc8/dfM/E4yKxbt06fffaZLrjgAhPKaWGio6VVq06e9fD7793PhMiZfRt3Zt/PPuPMvmbMXV3t3Fq4bh1n9jVr7toen/oatkLdzT332ZzZ98AB50GmnNnXnLnrew372Zl9bYZx6p++Z3beeefp73//e70nxfNHxcXFio6OVlFRkaKiorw2r8Ph0MqVK3XttdfKfupfBfAaemwu+ms+emwu+ms+X/a4sb+/PT5GZvHixXrwwQd16NChs6kPAADgrHm8a+nmm29WWVmZzjnnHIWHh9dJaMeOHfNacQAAAKfjcZBZvHixCWUAAAB4zuMgM3HiRDPqAAAA8JjHx8hI0v79+/XQQw9p3LhxKigokCT9+9//1u7du71aHAAAwOl4HGSysrLUp08fbdq0Se+8845KS0slSTt37tQjjzzi9QIBAAAa4nGQefDBB/XHP/5RGRkZCq49UY6kwYMHa+PGjV4tDgAA4HQ8DjK7du3SjTfeWGd5bGysjh496pWiAAAAGsPjIBMTE6Pc3Nw6y7dv366OHTt6pSgAAIDG8DjI3HLLLXrggQeUl5cnm82mmpoaff7555ozZ44mTJhgRo0AAAD18jjILFy4UD179lRiYqJKS0t13nnn6aqrrtLll1+uhx56yIwaAQAA6uXxeWSCg4P14osvat68edq1a5dKS0t14YUXqltDF+wCAAAwicdBplZiYqISExO9WQsAAIBHmnRCPAAAAH9AkAEAAJZFkAEAAJZFkAEAAJZ1VkGmT58+ysnJ8VYtAAAAHjmrIHPo0CE5HA5v1QIAAOARdi0BAADLOqsgc+WVVyosLMxbtQAAAHikySfEk6SVK1d6qw4AAACPsWsJAABYFkEGAABYFkEGAABYFkEGAABY1lkd7FtaWqqamhq3ZVFRUWdVEAAAQGN5vEXm4MGDuu666xQREaHo6Gi1bt1arVu3VkxMjFq3bm1GjQAAAPXyeIvMb37zGxmGob/97W+Ki4uTzWYzoy4AAIAz8jjI7Ny5U1u3blWPHj3MqAcAAKDRPN61dMkll3ChSAAA4Bc83iLz0ksvaerUqfr+++91/vnny263u63v27ev14oDAAA4HY+DzA8//KD9+/fr9ttvdy2z2WwyDEM2m03V1dVeLRAAAKAhHgeZO+64QxdeeKHefPNNDvYFAAA+5XGQ+e677/Svf/1L5557rhn1AAAANJrHB/sOHjxYO3fuNKMWAAAAj3i8Reb666/XrFmztGvXLvXp06fOwb6/+tWvvFYcAADA6XgcZKZOnSpJWrBgQZ11HOwLAACak8dB5ufXVgIAAPAVrn4NAAAsq1FbZJ555hlNmTJFoaGheuaZZ0479t577/VKYQAAAGfSqCCzaNEijR8/XqGhoVq0aFGD42w2G0EGAAA0m0btWjp48KDatm3r+n9DXwcOHPDowdPT03XJJZeoVatWio2N1Q033KA9e/a4jSkvL1daWpratm2ryMhIjR07Vvn5+R49DgAAaJk8PkZmwYIFKisrq7P8xIkT9X6S6XSysrKUlpamjRs3KiMjQw6HQ8OGDdPx48ddY2bNmqX3339fy5YtU1ZWlo4cOaIxY8Z4WjYAAGiBPA4y8+fPV2lpaZ3lZWVlmj9/vkdzrVq1SpMmTVLv3r3Vr18/vfLKK8rOztbWrVslSUVFRXr55Zf11FNPafDgwerfv7+WLl2q9evXa+PGjZ6WDgAAWhiPg0ztxSF/bufOnWrTps1ZFVNUVCRJrnm2bt0qh8Oh1NRU15iePXsqKSlJGzZsOKvHAgAA1tfo88i0bt1aNptNNptN3bt3dwsz1dXVKi0tdZ0srylqamo0c+ZMDRw4UOeff74kKS8vT8HBwYqJiXEbGxcXp7y8vHrnqaioUEVFhet2cXGxJMnhcMjhcDS5vp+rncubc8IdPTYX/TUfPTYX/TWfL3vc2MdsdJBZvHixDMPQHXfcofnz5ys6Otq1Ljg4WJ07d1ZKSornlf5HWlqavvrqK61bt67Jc0jOA4jr28W1evVqhYeHn9Xc9cnIyPD6nHBHj81Ff81Hj81Ff83nix7XdzxufRodZCZOnChJ6tKliy6//PI611g6G9OmTdMHH3ygtWvXqlOnTq7l8fHxqqysVGFhodtWmfz8fMXHx9c719y5czV79mzX7eLiYiUmJmrYsGGKioryWs0Oh0MZGRkaOnSoV3uBk+ixueiv+eixueiv+XzZ49o9Kmfi8SUKrr76alVXV2v58uX65ptvJEnnnXeeRo8eraAgz6YzDEPTp0/Xu+++q8zMTHXp0sVtff/+/WW327VmzRqNHTtWkrRnzx5lZ2c3uPUnJCREISEhdZbb7XZTvglmzYuT6LG56K/56LG56K/5fNHjxj6ex0Fm9+7d+tWvfqW8vDz16NFDkvSnP/1J7du31/vvv+86vqUx0tLS9MYbb+i9995Tq1atXMe9REdHKywsTNHR0Zo8ebJmz56tNm3aKCoqStOnT1dKSoouu+wyT0sHAAAtjMdB5s4771Tv3r21ZcsWtW7dWpL0008/adKkSZoyZYrWr1/f6Lmef/55SdKgQYPcli9dulSTJk2S5DyrcEBAgMaOHauKigoNHz5czz33nKdlAwCAFsjjILNjxw63ECM5P9H02GOP6ZJLLvFoLsMwzjgmNDRUS5Ys0ZIlSzwtFQAAtHAen0eme/fu9V4ioKCgQOeee65XigIAAGgMj4NMenq67r33Xi1fvlyHDx/W4cOHtXz5cs2cOVN/+tOfVFxc7PoCAAAwk8e7lkaNGiVJuummm1wnxavdRXT99de7bttsNlVXV3urTgAAgDo8DjKffvqpGXUAAAB4rEnnkQEAAPAHHgcZSSosLNTLL7/sOiFe7969dccdd7hdtgAAAMBsHh/su2XLFp1zzjlatGiRjh07pmPHjumpp57SOeeco23btplRIwAAQL083iIza9Ys/epXv9KLL77ouiRBVVWV7rzzTs2cOVNr1671epEAAAD18TjIbNmyxS3ESFJQUJB++9vf6uKLL/ZqcQAAAKfj8a6lqKgoZWdn11mek5OjVq1aeaUoAACAxvA4yNx8882aPHmy3n77beXk5CgnJ0dvvfWW7rzzTo0bN86MGgEAAOrl8a6l//3f/5XNZtOECRNUVVUlyXmp7bvvvluPP/641wsEAABoiMdBJjg4WE8//bTS09O1f/9+SdI555yj8PBwrxcHAABwOk06j4wkhYeHq0+fPt6sBQAAwCMeHyMDAADgLwgyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsnwaZNauXavrr79eCQkJstlsWrFihdt6wzA0b948dejQQWFhYUpNTdXevXt9UywAAPA7Pg0yx48fV79+/bRkyZJ61z/xxBN65pln9MILL2jTpk2KiIjQ8OHDVV5e3syVAgAAfxTkywcfOXKkRo4cWe86wzC0ePFiPfTQQxo9erQk6bXXXlNcXJxWrFihW265pTlLBQAAfsinQeZ0Dh48qLy8PKWmprqWRUdHa8CAAdqwYUODQaaiokIVFRWu28XFxZIkh8Mhh8Phtfpq5/LmnHBHj81Ff81Hj81Ff83nyx439jH9Nsjk5eVJkuLi4tyWx8XFudbVJz09XfPnz6+zfPXq1QoPD/dukZIyMjK8Pifc0WNz0V/z0WNz0V/z+aLHZWVljRrnt0GmqebOnavZs2e7bhcXFysxMVHDhg1TVFSU1x7H4XAoIyNDQ4cOld1u99q8OIkem4v+mo8em4v+ms+XPa7do3Imfhtk4uPjJUn5+fnq0KGDa3l+fr4uuOCCBu8XEhKikJCQOsvtdrsp3wSz5sVJ9Nhc9Nd89Nhc9Nd8vuhxYx/Pb88j06VLF8XHx2vNmjWuZcXFxdq0aZNSUlJ8WBkAAPAXPt0iU1paqn379rluHzx4UDt27FCbNm2UlJSkmTNn6o9//KO6deumLl266OGHH1ZCQoJuuOEG3xUNAAD8hk+DzJYtW3TNNde4btce2zJx4kS98sor+u1vf6vjx49rypQpKiws1BVXXKFVq1YpNDTUVyUDAAA/4tMgM2jQIBmG0eB6m82mBQsWaMGCBc1YFQAAsAq/PUYGAADgTAgyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsggyAADAsiwRZJYsWaLOnTsrNDRUAwYM0BdffOHrkgAAgB/w+yDz9ttva/bs2XrkkUe0bds29evXT8OHD1dBQYGvSwMAAD4W5OsCzuSpp57SXXfdpdtvv12S9MILL+jDDz/U3/72Nz344IONnud45XEFVgbWWR4YEKjQoFC3cQ0JsAUozB7mul1eXa7jlcdlN+xnHFvmKJNhGPXOa7PZFG4Pb9LYE44TqjFqGqw5IjiiSWPLq8pVXVPtlbHh9nDZbDZJUkVVhapqqho99nQ9DrOHKcDmzOKV1ZVyVDsanNeTsaFBoQoMCPR4rKPaocrqygbHhgSFKCggyOOxVTVVqqiqaHBscGCw7IF2j8dW11Sftr/2QLuCA4NPjq0qb3DeU8fWGDU64TjhlbFBAUEKCQqRJBmGoTJHmVfGevJzfzbvEccrjzfYY94jTuI9wskf3yNO9xpujveIxrAZDf00+IHKykqFh4dr+fLluuGGG1zLJ06cqMLCQr333nt17lNRUaGKipPfpOLiYiUmJkoPSgqtM1wjzxmp924+OU/M/8Q0+AZ4VdJV+vg3H0uSHA6HOjzVQcXVxfWO7d+hvzbcvsF1u9uSbvqu6Lt6x/Zq10s7p+x03e7313765ug39Y5Njk7W3rS9rtspS1O0NXdrvWPbhbXTkVlHXLdT/y9Va7PX1js23B6uwvsLXbdHvz1a/97/73rHSlLl707+gN3yzi1659t3Ghz705yfXG9qk9+frNd3vd7g2O9nfK/2Ee0lSdNWTtNfd/y1wbH/757/p84xnSVJD655UE9teqrBsdvv2q7e7XtLkhasXaA/rvtjg2PXT1qvixMuliQ9ufFJzf1kboNjM8Zn6OrkqyVJz295XjNWz2hw7IqbVujac6+VJL325Wu684M7Gxz7xo1v6L96/Zckafk3y3Xru7c2OPalUS9pQt8JkqSV+1bqhn/c0ODYp4c9rbsvvluStGb/Go18e2SDY9MHp+u+y+6TJG05skWXv3J5g2MfuuIhzbtqniRp9w+7deGLFzY4dvaA2Xp8yOOSpEOFh9T9ue4Njp160VQ9M+IZSdIPx39Qx6c7Njj2tj636eXrX5bkDBCt/7d1g2PH9Byjt8a85bodvLDhN82mvkdIUsKiBB09cbTesbxHnMR7hJO/vUdkfZeloX8f2uBYs98jiouL1a5dOxUVFSkqKqrB8X69Rebo0aOqrq5WXFyc2/K4uDh9++239d4nPT1d8+fPb/RjFPxQoJUrV7puV1c3/FfDjz/+6Db2dIoKi9zGlpU1/NdhaWmp29jS0tIGx5aVlbmNLSosanBsZWWl29gff/yxwbHV1dVuYwt+OP2uu1PH5uXmnXbsRx99pNBAZ4o8fPjwacd+vOZjRQdFS5JyDuecduynn36quBDna+PA9wdOO/aztZ/puzDnL4m9uXtPO/bz9Z+rINz5/L8tqP91Vmvjxo06vtv5F/ruH3afduyWzVuk/+f8/84fd5527Pbt2xV+0PlX9fbC7acdu3PnTq087Px+bCnactqxu3fv1soC59hdJbtOO/bbb7/VymPOsXvLTt+zvXv3amWpc2z2iezTjj1w4IBWVjjH5lfkn3bsd9nfuV5rRVUNv9Yl52urdmx5dcN/GUrO12xjf5bP5j2isrLhv6h5jziJ9wgn3iOcat8jTvczcSq/3iJz5MgRdezYUevXr1dKSopr+W9/+1tlZWVp06ZNde7T0BaZ7458V2+ia+pmY4fDofdXva/BgwfLbmfXkhmbjUtPlGr1x6sb7DGbjZ2autm4vKJcK1evbLC/7FryfOzPf+4Ljxfqk08+qbfHvEecxHuEk7+9R1TXVKvkREmDr2Gz3yNaxBaZdu3aKTAwUPn57n+15efnKz4+vt77hISEKCQkpM7ymIgYRUU03AjXOHtMo+sLDQxVTERMvT9APxdtj270vJ6MbcxjW3VspCIb3WN/qNdutytc4Wce6OlY2RUWEnbmgR6OlRr/GrbLrtCQevbNNiAkuO7PoDfGBgc3fr+5J2M9+bn3aGxETKN7zHuE52N5j/jPWJPeI+yyO4N8I3psxntEY3vr159aCg4OVv/+/bVmzRrXspqaGq1Zs8ZtCw0AAPhl8ustMpI0e/ZsTZw4URdffLEuvfRSLV68WMePH3d9igkAAPxy+X2Qufnmm/XDDz9o3rx5ysvL0wUXXKBVq1bVOQAYAAD88vh9kJGkadOmadq0ab4uAwAA+BlLBJmzUXtkf3Fx/ed7aSqHw6GysjIVFxd7dLAXGo8em4v+mo8em4v+ms+XPa79vX2mD1e3+CBTUlIiSc6T4gEAAEspKSlRdHTDn9Tz6/PIeENNTY2OHDmiVq1auc494A2156fJyck57efb0XT02Fz013z02Fz013y+7LFhGCopKVFCQoICAhr+kHWL3yITEBCgTp06mTZ/VFQUP0Amo8fmor/mo8fmor/m81WPT7clppZfn0cGAADgdAgyAADAsggyTRQSEqJHHnmk3sshwDvosbnor/nosbnor/ms0OMWf7AvAABoudgiAwAALIsgAwAALIsgAwAALIsgAwAALIsg00RLlixR586dFRoaqgEDBuiLL77wdUmW9Oijj8pms7l99ezZ07W+vLxcaWlpatu2rSIjIzV27Fjl5+f7sGL/t3btWl1//fVKSEiQzWbTihUr3NYbhqF58+apQ4cOCgsLU2pqqvbu3es25tixYxo/fryioqIUExOjyZMnq7S0tBmfhf86U38nTZpU5zU9YsQItzH0t2Hp6em65JJL1KpVK8XGxuqGG27Qnj173MY05n0hOztb1113ncLDwxUbG6v7779fVVVVzflU/FJj+jto0KA6r+GpU6e6jfGn/hJkmuDtt9/W7Nmz9cgjj2jbtm3q16+fhg8froKCAl+XZkm9e/dWbm6u62vdunWudbNmzdL777+vZcuWKSsrS0eOHNGYMWN8WK3/O378uPr166clS5bUu/6JJ57QM888oxdeeEGbNm1SRESEhg8frvLycteY8ePHa/fu3crIyNAHH3ygtWvXasqUKc31FPzamforSSNGjHB7Tb/55ptu6+lvw7KyspSWlqaNGzcqIyNDDodDw4YN0/Hjx11jzvS+UF1dreuuu06VlZVav369Xn31Vb3yyiuaN2+eL56SX2lMfyXprrvucnsNP/HEE651ftdfAx679NJLjbS0NNft6upqIyEhwUhPT/dhVdb0yCOPGP369at3XWFhoWG3241ly5a5ln3zzTeGJGPDhg3NVKG1STLeffdd1+2amhojPj7e+J//+R/XssLCQiMkJMR48803DcMwjK+//tqQZGzevNk15t///rdhs9mM77//vtlqt4Kf99cwDGPixInG6NGjG7wP/fVMQUGBIcnIysoyDKNx7wsrV640AgICjLy8PNeY559/3oiKijIqKiqa9wn4uZ/31zAM4+qrrzZmzJjR4H38rb9skfFQZWWltm7dqtTUVNeygIAApaamasOGDT6szLr27t2rhIQEde3aVePHj1d2drYkaevWrXI4HG697tmzp5KSkuh1Ex08eFB5eXluPY2OjtaAAQNcPd2wYYNiYmJ08cUXu8akpqYqICBAmzZtavaarSgzM1OxsbHq0aOH7r77bv3444+udfTXM0VFRZKkNm3aSGrc+8KGDRvUp08fxcXFucYMHz5cxcXF2r17dzNW7/9+3t9af//739WuXTudf/75mjt3rsrKylzr/K2/Lf6ikd529OhRVVdXu30DJSkuLk7ffvutj6qyrgEDBuiVV15Rjx49lJubq/nz5+vKK6/UV199pby8PAUHBysmJsbtPnFxccrLy/NNwRZX27f6Xr+16/Ly8hQbG+u2PigoSG3atKHvjTBixAiNGTNGXbp00f79+/W73/1OI0eO1IYNGxQYGEh/PVBTU6OZM2dq4MCBOv/88yWpUe8LeXl59b7Ga9fBqb7+StKtt96q5ORkJSQk6Msvv9QDDzygPXv26J133pHkf/0lyMCnRo4c6fp/3759NWDAACUnJ+sf//iHwsLCfFgZ0DS33HKL6/99+vRR3759dc455ygzM1NDhgzxYWXWk5aWpq+++srtuDl4T0P9PfV4rT59+qhDhw4aMmSI9u/fr3POOae5yzwjdi15qF27dgoMDKxzhHx+fr7i4+N9VFXLERMTo+7du2vfvn2Kj49XZWWlCgsL3cbQ66ar7dvpXr/x8fF1DlyvqqrSsWPH6HsTdO3aVe3atdO+ffsk0d/GmjZtmj744AN9+umn6tSpk2t5Y94X4uPj632N165Dw/2tz4ABAyTJ7TXsT/0lyHgoODhY/fv315o1a1zLampqtGbNGqWkpPiwspahtLRU+/fvV4cOHdS/f3/Z7Xa3Xu/Zs0fZ2dn0uom6dOmi+Ph4t54WFxdr06ZNrp6mpKSosLBQW7dudY355JNPVFNT43pDQ+MdPnxYP/74ozp06CCJ/p6JYRiaNm2a3n33XX3yySfq0qWL2/rGvC+kpKRo165dboExIyNDUVFROu+885rnifipM/W3Pjt27JAkt9ewX/W32Q8vbgHeeustIyQkxHjllVeMr7/+2pgyZYoRExPjdgQ3Gue+++4zMjMzjYMHDxqff/65kZqaarRr184oKCgwDMMwpk6daiQlJRmffPKJsWXLFiMlJcVISUnxcdX+raSkxNi+fbuxfft2Q5Lx1FNPGdu3bze+++47wzAM4/HHHzdiYmKM9957z/jyyy+N0aNHG126dDFOnDjhmmPEiBHGhRdeaGzatMlYt26d0a1bN2PcuHG+ekp+5XT9LSkpMebMmWNs2LDBOHjwoPHxxx8bF110kdGtWzejvLzcNQf9bdjdd99tREdHG5mZmUZubq7rq6yszDXmTO8LVVVVxvnnn28MGzbM2LFjh7Fq1Sqjffv2xty5c33xlPzKmfq7b98+Y8GCBcaWLVuMgwcPGu+9957RtWtX46qrrnLN4W/9Jcg00bPPPmskJSUZwcHBxqWXXmps3LjR1yVZ0s0332x06NDBCA4ONjp27GjcfPPNxr59+1zrT5w4Ydxzzz1G69atjfDwcOPGG280cnNzfVix//v0008NSXW+Jk6caBiG8yPYDz/8sBEXF2eEhIQYQ4YMMfbs2eM2x48//miMGzfOiIyMNKKioozbb7/dKCkp8cGz8T+n629ZWZkxbNgwo3379obdbjeSk5ONu+66q84fOfS3YfX1VpKxdOlS15jGvC8cOnTIGDlypBEWFma0a9fOuO+++wyHw9HMz8b/nKm/2dnZxlVXXWW0adPGCAkJMc4991zj/vvvN4qKitzm8af+2gzDMJpv+w8AAID3cIwMAACwLIIMAACwLIIMAACwLIIMAACwLIIMAACwLIIMAACwLIIMAACwLIIMAACwLIIMAL83adIk3XDDDb4uA4AfIsgAAADLIsgA8BvLly9Xnz59FBYWprZt2yo1NVX333+/Xn31Vb333nuy2Wyy2WzKzMyUJOXk5Oimm25STEyM2rRpo9GjR+vQoUOu+Wq35MyfP1/t27dXVFSUpk6dqsrKSt88QQBeF+TrAgBAknJzczVu3Dg98cQTuvHGG1VSUqLPPvtMEyZMUHZ2toqLi7V06VJJUps2beRwODR8+HClpKTos88+U1BQkP74xz9qxIgR+vLLLxUcHCxJWrNmjUJDQ5WZmalDhw7p9ttvV9u2bfXYY4/58ukC8BKCDAC/kJubq6qqKo0ZM0bJycmSpD59+kiSwsLCVFFRofj4eNf4//u//1NNTY1eeukl2Ww2SdLSpUsVExOjzMxMDRs2TJIUHBysv/3tbwoPD1fv3r21YMEC3X///frDH/6ggAA2SgNWx08xAL/Qr18/DRkyRH369NGvf/1rvfjii/rpp58aHL9z507t27dPrVq1UmRkpCIjI9WmTRuVl5dr//79bvOGh4e7bqekpKi0tFQ5OTmmPh8AzYMtMgD8QmBgoDIyMrR+/XqtXr1azz77rH7/+99r06ZN9Y4vLS1V//799fe//73Ouvbt25tdLgA/QZAB4DdsNpsGDhyogQMHat68eUpOTta7776r4OBgVVdXu4296KKL9Pbbbys2NlZRUVENzrlz506dOHFCYWFhkqSNGzcqMjJSiYmJpj4XAM2DXUsA/MKmTZu0cOFCbdmyRdnZ2XrnnXf0ww8/qFevXurcubO+/PJL7dmzR0ePHpXD4dD48ePVrl07jR49Wp999pkOHjyozMxM3XvvvTp8+LBr3srKSk2ePFlff/21Vq5cqUceeUTTpk3j+BighWCLDAC/EBUVpbVr12rx4sUqLi5WcnKynnzySY0cOVIXX3yxMjMzdfHFF6u0tFSffvqpBg0apLVr1+qBBx7QmDFjVFJSoo4dO2rIkCFuW2iGDBmibt266aqrrlJFRYXGjRunRx991HdPFIBX2QzDMHxdBACYYdKkSSosLNSKFSt8XQoAk7BtFQAAWBZBBgAAWBa7lgAAgGWxRQYAAFgWQQYAAFgWQQYAAFgWQQYAAFgWQQYAAFgWQQYAAFgWQQYAAFgWQQYAAFgWQQYAAFjW/wddyqtK/KuhRgAAAABJRU5ErkJggg==\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "file_names = []\n",
        "for i in tqdm(range(T_all)):\n",
        "  file_name = f\"fig_data/progress_{i:04d}.png\"\n",
        "  file_names.append(file_name)\n",
        "  plot_progress(y_best_data[:i+1], y_data[init_data_num-1:init_data_num+i], optim, file_name)"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "064nJwu31zoA"
      },
      "source": [
        "すると，あとは以下のように実行するだけでアニメーションとして見ることができます．  \n",
        "縦軸が最適解までとのエネルギーの差です．つまり，縦軸が0だと最適解を見つけることができていることになります．  \n",
        "横軸はステップ数です．先ほど記述した一連のアルゴリズムを一回繰り返すことが1ステップに相当します．つまり1ステップにつき線形回帰するためのデータセットが1つ追加されるという感じです．  \n",
        "赤のバツ印がそのステップの時点における代理モデルの最適解です．  \n",
        "そして，青い実線が現在持っているデータセットの中で最小のエネルギーを示しています．  "
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 488
        },
        "id": "-Hui3l-UiNEL",
        "outputId": "00d3de10-d10d-4a54-d093-2b1eab5caccb"
      },
      "outputs": [],
      "source": [
        "APNG.from_files(file_names, delay=100).save(\"animation.png\")\n",
        "Image(\"animation.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "metadata": {
        "id": "nAsup4oI4z6g"
      },
      "source": [
        "さて，これを見てもらうといかがでしょうか．もちろんブラックボックス関数のQUBOによって異なりますが，早い段階で収束してるのではないでしょうか．これはある時点から同じデータしかサンプリングしかせず，データ数が増えても係数を推定するのに必要な情報は増えておらず，代理モデルが更新されない現象が起きているためです．  "
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "### 後処理の追加"
      ],
      "metadata": {
        "id": "5EoU5z7rVY9q"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "データセットに重複したものがサンプリングされるのが良くないのでしたら，重複した時は代わりにランダムに生成したものを$\\vec{x}^{next}$として$y=f(\\vec{x}^{next})$を計算し，データセットに追加するようにプログラムを少し書き換えてみましょう．"
      ],
      "metadata": {
        "id": "f7eKx1GsVVn4"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "id": "8RoALdMnMCuZ"
      },
      "outputs": [],
      "source": [
        "# x_dataとy_dataをリストに変換\n",
        "x_data = init_x_data.tolist()\n",
        "y_data = init_y_data.tolist()\n",
        "\n",
        "# zを計算しておく\n",
        "z_data = []\n",
        "for i in range(len(x_data)):\n",
        "  z = calc_z_from_x(x_data[i])\n",
        "  z_data.append(z)\n",
        "\n",
        "# yの中で最もエネルギーが低い値を保存\n",
        "y_data.sort(reverse=True)\n",
        "y_min = y_data[-1]\n",
        "y_best_data = []\n",
        "y_best_data.append(y_min)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "T1Y3I_WqmqfF",
        "outputId": "ec2a0fe1-a3b7-40fd-8a4b-c667dc51eaad"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "100%|██████████| 256/256 [00:37<00:00,  6.88it/s]\n"
          ]
        }
      ],
      "source": [
        "for t in tqdm(range(T_all)):\n",
        "  # 1. 代理モデルを作成する\n",
        "  # データセットから最もフィットする係数ベクトルaを求める\n",
        "  a = predict_a(z_data, y_data, alpha=0.01)\n",
        "  # 係数ベクトルaからQUBOに変換する\n",
        "  Q = calc_Q_from_a(a)\n",
        "\n",
        "  # 2. QUBOを最適解を求める\n",
        "  x_next = sample_from_qubo(sampler, Q)\n",
        "\n",
        "  # この部分を先ほどのプログラムに書き加える\n",
        "  # x_nextがすでにデータセットに含まれている場合はランダムにサンプリング\n",
        "  if any([x_next == x for x in x_data]):\n",
        "    x_next = np.random.randint(0, 2, N).tolist()\n",
        "\n",
        "  # 3. y_next=f(x_next)を求める\n",
        "  y_next = black_box(x_next)\n",
        "\n",
        "  # 4. データセットに追加する\n",
        "  x_data.append(x_next)\n",
        "  y_data.append(y_next)\n",
        "\n",
        "  z_next = calc_z_from_x(x_next)\n",
        "  z_data.append(z_next)\n",
        "\n",
        "  # 現時点のデータセットで最もエネルギーが低いyをy_best_dataに追加\n",
        "  if y_next < y_best_data[-1]: # 一番最後尾の値と比較\n",
        "    y_best_data.append(y_next)\n",
        "  else:\n",
        "    y_best_data.append(y_best_data[-1])"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/"
        },
        "id": "Ga6Udi22LsDW",
        "outputId": "39805b98-a973-4f84-e803-daf08fe840ef"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stdout",
          "text": [
            "mkdir: cannot create directory ‘fig_data_random’: File exists\n"
          ]
        }
      ],
      "source": [
        "!mkdir fig_data_random"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "さて，結果をプロットしてみましょう．"
      ],
      "metadata": {
        "id": "noBJDNOrb73g"
      }
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 472
        },
        "id": "Ed-DmxDLGN3r",
        "outputId": "317acbb0-7fc0-493b-b987-f0cfc54f8110"
      },
      "outputs": [
        {
          "output_type": "stream",
          "name": "stderr",
          "text": [
            "100%|██████████| 256/256 [00:52<00:00,  4.85it/s]\n"
          ]
        },
        {
          "output_type": "display_data",
          "data": {
            "text/plain": [
              "<Figure size 640x480 with 1 Axes>"
            ],
            "image/png": "iVBORw0KGgoAAAANSUhEUgAAAjIAAAG2CAYAAABlBWwKAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjcuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/bCgiHAAAACXBIWXMAAA9hAAAPYQGoP6dpAACGxklEQVR4nO2deXxU1fn/P5M9ISRhzUIyhMgqAlpQjFAWRTYXFH5VARFUtChuULvQWgW7YG2r1hb1W7VqaUGrRSw2gogSZFVQdklZEjIhCQExG4FkkpzfHyd35s7MnZl7Z+4687xfr7wmc++ZO+eeOefcz3nOc55jY4wxEARBEARBWJAYozNAEARBEAQRKiRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLIYKmaVLl8Jms3n8DRw40HX+4sWLWLhwIbp164bU1FTMmDEDp0+fNjDHBEEQBEGYCcMtMoMHD0ZVVZXrb+vWra5zixYtwrp16/Duu++iuLgYlZWVmD59uoG5JQiCIAjCTMQZnoG4OGRlZfkcr6urw+uvv45Vq1bh2muvBQC88cYbGDRoEHbu3Imrr75a76wSBEEQBGEyDBcyR48eRU5ODpKSklBYWIjly5fDbrdjz549cDqdmDBhgivtwIEDYbfbsWPHDr9Cprm5Gc3Nza737e3tOHfuHLp16wabzab5/RAEQRAEET6MMTQ0NCAnJwcxMf4nkAwVMiNHjsSbb76JAQMGoKqqCsuWLcP3v/99HDx4ENXV1UhISEBGRobHZzIzM1FdXe33msuXL8eyZcs0zjlBEARBEHrgcDiQm5vr97yhQmbKlCmu/4cOHYqRI0eid+/e+Ne//oXk5OSQrrlkyRIsXrzY9b6urg52ux2lpaXo3Llz2HkWcDqd+OyzzzB+/HjEx8erdl3daWgAGhuB7Gzfc1VVQGoqoGK5KcGyZWziMhWjuHwrK4H/9/+A8nL3MbsdeO89ICdHu4xaGMkyDlaOVVXAjBn8vPi4+HN2O/Dvf0vXMa3QIl9K65RXemdyMj5bsQLjBw1CfF6eCjcZAg0NwOzZwNmzvvkW8tu9O/DPf5qi3SvFyH64oaEBffr0CfrsNnxqSUxGRgb69++PY8eO4frrr0dLSwtqa2s9rDKnT5+W9KkRSExMRGJios/xrl27Ii0tTbW8Op1OpKSkoFu3btZ6yIqpqwNmzgRqaoDNmwFxR+BwANOmAT17AuvXA+npumfPsmXcrVto53RGUfk6HMAttwAnTgAFBcDKlcCcOcD//sePe9cfAoCfMu7WDXjjDWDUKHfCN94Ahgzh/8fFAV27Aq2twCefuMu1Wzf+ftw4fj4vT3m7rKvjD16p0W1FBX/Q+rumFvkKVhZiJOqgc/58Xr6zZiF+wwZj6mC3bkBRkXS5dusGbNoUuFxNjpH9sPB9wdxCDF+1JKaxsRHHjx9HdnY2hg8fjvj4eGzatMl1vqSkBOXl5SgsLDQwlxFEQwMXMSdO8E7I4eDHHQ7+/sQJfr6hwchcEkZTUeGuDwUFXLRccw1/LShw15+KCmPzaRUcDi4CxcyZ425/6el88FBc7Ptgzsvjx0MZXNTVAZMnA2PHur9LnKexY/n5ujrpz2uRr2BlIeCvDhYV8fNlZcbWwfR0aXEI8OMWFTFWwVAh8/jjj6O4uBhlZWXYvn07br31VsTGxmLmzJlIT0/Hvffei8WLF+Ozzz7Dnj17cPfdd6OwsJBWLKlFbq7vw2j7dt8OI8DcJBEFdO7MLXNCfRAeYnl57vrTs6clzea6Ix4kFBQA27Z5tj+xmFH7wajGwCXUfNXV+YoM8ffm5/svC8B/HezVi7/m51MdjGIMnVqqqKjAzJkz8e2336JHjx4YPXo0du7ciR49egAAnn/+ecTExGDGjBlobm7GpEmT8NJLLxmZ5chDeBgJHYpg4vXuMIjoRRiJS5nOhZG4hU3nuiFlVfBuf+PG8fLUYvAgDFzE3yVMEWo5cBEsQeIpbHFZxMUBXboAgwf7L4tAdRAAPvoIyMigOhilGCpk3n777YDnk5KSsGLFCqxYsUKnHEUpeXm8QxPPU69cSSKGcJOe7v8hEe6DLxy/DS1RO1+CVQGQtmyNG6e9VcGIgYu3JWjzZiAtjQuPuDjucyOUdaCyCFQHc3IAhf4bbW1tcDqdod2T1WloAM6fB6T8TaurgU6dXGXvdDoRFxeHixcvoq2tTdVsxMfHIzY2NuzrmMrZlzAIf/PUkWqRMeuDMxqRGq0LCFMPRjica5Evs1i2Qhm4hNNm/FmCzp3jIsbbEqRxWTDGUF1djdraWtWvbQna23m9bmvjv2mcSAa0tgKnTwOxsbx+x8SAMYasrCw4HA5NYrFlZGQgKysrrGuTkIl2vOfsxaZmYfQUSWLGrA/OaEVqtJ6X51kvhXR6/h5a5UtLy5ZclA5c1GgzSi1BGpaFIGJ69uyJlJSU6AuU2tLCRUxLC8AY9zNKSODvy8q4pSwhgf8mCQlob29HY2MjUlNTAwalUwpjDE1NTaipqQEAZIcRSoCETDRj9Jy9XOrrgQsX1LGgmPXBaWW0GK1r7bcRDLPmK1xCGbio1WZMMIXd1tbmEjHdTBQKQVeSkoCBA4GSEqC5GTh5EujTh7+2tACJicCAAVzMgEfHb2lpQVJSkqpCBoArXlxNTQ169uwZ8jSTqZZfRw1SHvwCFRX+lz+qjVVWo0yfHvqSUW9opZa6hLukF/Csb8Jo3VtcG4FZ8xUqoS6jV6vNyF1qrSGCT0xKSoq8D7S28oe7FC0t/LwVSUjgYiUxkYuZI0f4q5eI0QPhtwjHX4mEjN6o0fGrhVbxKtTmzBl1Y91E2gPKSNSKRSSM1sWYweHcX77S0swxGJFLXR23bEoNXGw24D//CTxwCbfNyF12rhOyppNaW4GjR7nlwlvMtLTw40ePWlvM9OnjeaxPH11FDCDztwgCCRm9MVsQOisEcioqUt+CYtYHp9XwHq2PHi39+wCBH+4mGK1LIpWv2bOBa681x2BEDsLg6aabgNde8xy4CPmdPx9Yty64n0sobcaqARXb27lIaW72FDOCiGluBpxOns4bK1hrWlqA0lLPY6Wl/i1QJoaEjN7Q1IZyevVS34Ji1genFcnL4yP6xES+B4737wMEfribbLQeNF9lZcD+/eYZjARDPHi6+Wbu4An45jctLfDAJdQ2Y5UpbG+8p19KSvgeaoKI8WdJCMdaE2Aqa9zYsXjskUcU3oQfxGIsMZH7zIjv02JihoSMEdDUhnLUtKCY9cFpZdLTffeREn6vQA93s47Wg+WrtZUvW7XCYESNwVM4bcYqU9hS+PMlSUjgcWsEQSBlrWltlbbW+CPYVFZTE1BbG76lx1vEDBiAzbt3wzZ0KGqbmy0pZkjIGAVNbShDLQuKWR+cVocxz3gUAPCDH/CpJnFZd+7sWbbi0frq1dwqABg/WpdjRRgyhIfGt8JgJJzBkxptRu4UtlkWQoiR8iUpKOBWjIQEaWtNYiJPo2SVT7CprPZ29184xMTwtirl2NuvHz8eF6cs7wZjnZxGGjS1IZ9Tp9SzoFjVzG1mhNF6eTlgt/MoqwBQWek+tnkzFyneju7CaH31ar4Tu3j6ycjRupCvdet8pxCEfP3jH8D//Z/nOTMPRkIdPOnVZsy0EEKMP1+Stjb+wLfZfFf+9OnD+6f//Y+HjvB3XbF1JdhUVkwMWjt1wkOPPYb09HR0794dv/zlL8EYA1pa0Hz+PB5ftAi9evVCp06dMHLkSGwWpnYBnDx5EjfddBO69OiBTiNGYPDtt6Pok09QVlaG8ePHAwC6ZGbCNnQo5i1f7jswMTHWyWkkEW1B6MJl6lT1Yt2YJbpqpCA1Wnc4PGOFtLbyzt5fPJL6ei5ipOKRGD09c++90oHgGAOmTOFRUMWYOSJ2qBG89WozBsV4YozP2kjS0sLFSIsTSEjiFriyMqC2Bdh3jKdx2gCIxG5WLnCojH/G1grUH/O1fHRcN6VzLGz9+7lFgyBmBPFy5Ij7uM2Gt/7xD9w7fTq+WLkSu7/7Dvc/+CDsOTm4b/RoPLRsGQ4fP463n3kGOVddhffXrcPkyZNx4MAB9OvXDwsXLkRLSwu2bNmCTp064fDhw0hNTUVeXh7+/e9/Y8aMGSgpKUFaWporvotlYBFOXV0dA8Dq6upUvW5LSwtbu3Yta2lpUfZBh4OxggLGAP5aXs6Pl5d7Hnc4VM2vFXGV8ZgxnmUlIJTZ1VczVltrTCYtTMh1WExtLS9/4fcR12PhLzGRsYMHeXrver5tm3R7MAOB2qrd7r4/u93vfahSxmpglXJXmE+l5XvhwgV2+PBhduHCBdexxkbP6qrnX+POA4w1N/tmtKGBsS+/dP+dO8fGjhjBBvXpw9qFY/v3s58uWsQGFRSwk+vWsdjYWHaqqIix/ftd17zuuuvYkiVLGGOMDRkyhC1dulSyXD777DMGgH333Xc+59ra2th3333H2traZJWxUqR+EwG5z2+aWtIbmtpQzpo11nQUjAbEjpyA7xSg3c5HljffzEfWVnJ09+ckO3o0nzID+P1t3WpuPysr+YVZqX6oQf/+vnFbpKayKiqApCRcPWwYt/10TGcV9uqFoydP4sDx42hra0P/GTOQWliI1K5dkZqaiuLiYhw/fhwA8Mgjj+DXv/41Ro0ahaeeegr79+/X5Rb1gKaW9IamNpSTlua7IkbA6KkHgtfVhgbp7S62bvWdAtQ6VL2am4L62yMoMRHIzOT3Z8Qu1koww67bStB5K4OUFO6KIomw8khKbJSU8DgyCQk8RIS3+EhI4EIF6JieavE5l5Lh57piX5vSUv7+4kUgNdXtQyOi8fx5xMbGYs+uXYj1mhZKTU0FAMyfPx+TJk3Cf//7X3z88cdYvnw5/vjHP+Lhhx+WU0ymhoSMEZhh4ziCUBMlD0std1vXYlNQqQfrf/4DXHpp8MFIGGHXVcNqgyct64cENhvQqZO/s34ekYkxQHoc0AqgTz4XG8mi1UQ2G5AcA3TuWCE0ON/t7wIAA/OB1CAiRvCrEXxm2tuxa8cO/ht2WFl2HjiAfnY7rhgwAG1tbaiprcX3hwzxe695eXlYsGABFixYgCVLluDVV1/Fww8/jIQOodbW1hagpMwLTS0RBBE+cmOF1NdrG8NHi8jZUg/WBx5wB5bzxqiI2IGWLjc0+Le4mCWCN2CdGE9xcXypckGB22IiDizHmNuaIzeCrr9l0YKYiYlBeVUVFj/yCErKyrB6wwb8+V//wqN33IH+vXtj9g034K45c7BmzRqUlpbiiy++wPLly/Hf//4XAPDYY49hw4YNKC0txVdffYXPPvsMgwYNAgD07t0bNpsNH374Ic6cOYNGvyYqc0JCRg3MGPuAiF70rI/i7/KOFSL+rtxc6ekntX011I6cbZUHqxpLl43ux7Tw5Ql0T+FuI9DezvMktqCkprqXULe0cCuMeFl2oAi6gjiS2rQxIQFISsJdN92ECxcu4Kp587Dw2Wfx6B134P5bbwVsNrzxxBO4a/Jk/GjxYgwYMAC33HILvvzyS9jtdgDc2rJw4UIMGjQIkydPRv/+/fHSSy8BAHr16oVly5bhZz/7GTIzM/HQQw+FXi5GoIkbsonQfNXS2bOeqzbEqLGqprbW/womhyOiV+uYZsWHWujxW3qvIhLjVR/DLl8F3xVS+nCQWj2ldHWOCisMdavD4eZVz9/GHyHkwad8xW3M+3oOB7tQU8NXyNTV8dU9hw8z5nSGll+nk39etErIRXMzY/v2MbZ7t2uFkStNczN/7308EOLPiK/Z0CB9XM41ZUKrlqKBxkbtNoE0a4AoQjl6/ZZ6bkqq9LuCTT+tWwe8/rr0NIdSi4AakbOttMIwXEuUGTazDXcrA+82Jr6n0aO5j9N993G/pbKy0LYREBPMgtKvH/ck9jdVpCSCrjDtlJAAJCdLW4CSk/l5i0XlVQVNJJaJ0CWOjFYxGqI85oxqo1kzWLX0/C1l1kdVyletuq+2RUANi4yQrzDqju5WxXDu28yxZvz8Dh6Wcak25hXz58KVV7LDH3/MLmhguZDE6fT/Hc3NyqxBwrWkrik+HqqFyQ9WsMiQkAkRnw5KrY7TGzN3LhqjScA2MXoH1NPzt5RRH2WXb7CH+cGDppjGcWGiNmPI9Oi2bZ6/xbZt8j+rVT8WDgHacEtpqTtoZm2t9G8vFjK9e7PDH33ELuzbp72IiRBIyJgAXSP7htOBBMKMnYsOqPIQMJtVS8/fMkh9lFW+coXghg3h1301BIjJfm9LWWQEtOrHQiXAb9oyaBAv30GD3L+pVBnk5HgKmbNnjbsfi2EFIRNlE2kaouUmkLRTduiovYolXPT6LdWqj3J8Jyorue9BuN+lRlRXK/m1qI0aK6zMuJltoDZcVsbTFBW527BUG/Pm1Cnf5c+EddFEYpkIS/vICIQ7yjKDj0gIqDqaNYtVS498qO0jE+h6drvbdK9W3Q/XImCi+m6ZVUtSac02jS3RdlwWGXH5SrWxjrp6Yft2dnjjRv18ZCIAsshEA5WV2sbGCHeURSufOGawaukRk0SLWBz+LCUd8SlQXq7ed6lhEfCOZyPGTMHf1CRcS5QV9mOSasN//avne3Ebs9uBnBzP81lZfGuJhATpWC6EJSEhEy6pqdqZstXoXMywrNIMGG0y1+tBodXUitRD5NVX+YNCre+ySvA5MxLu0mWh3uTnA6tXS/+W6emhL1VWA6k2fP/97v+921hRERczwl95OXDXXTxtfr6y5c+EudHEVmQidJla0sqUHczRMj+fsSuukL6++HvNbjL2g2pmeTPcv56rp2TWR0Xl629K7OBBdeq+yZx01cJSQR1PnmTse9+TrqO7dvH+Rq8Vft74acMtycm8fEtLpduY0BY6Pn9h2jR2+OBBPo2hwVLlSISmlqIFrUzZgUZZaWlA167Ad9/x/WvEeE8ZqeFEaVXMYjIPd8Ss9LvUrI+BLCU338xlTbjfFUlOukaH+g+VmBigtlbaejtzJnesNcJ6G6gN5+fzNFOn8nx5tzGhLQht7NVX3RYYIXgcERbz5s3DLbfcYmgeSMiYHX8PpYYG3umUlcmbMjKDj4gRmOkBaUXfDb2EoJ5CT0us7JNmthV+AoHacFER/79HD/cO3oHamBWEMKEYEjJWRWmnY7SPiFFEygNSDlpYAvQUglYUet54+6SdOsWPnzplDZ80M1pvA7XhXr3465o1+tUPs1jcWlv9OyqHuyGmxSAhY2XkdjrR7kQZCQ/IYEhZAoQOV8oSUFHhOyUpRTQJQTXwHmBMncqPT51qrFVDCWpYb9V+2AdqwwCfatcDgy1u7733HoYMGYLk5GR0694dE8aMwfnvvsOXX36J66+/Ht27d0d6ejrGXnMNvvrgAw8xY7PZ8H//93+48cYbkZKSgkGDBmHHjh04duwYxo0bh06dOuGaa67B8ePHXZ9ZunQpvve97+GNN95A7969kZKSgttuuw11Ae6vvb0dy5cvR58+fZCcnIxhw4bhvffe06Q8BEjIqIkRSj1Yp2MWHxFCW7wtAYcO8Q511Ci+YZ7YEiB0uNOny7t2NAhBNREPMISAbWVl1vFJC9d6a+XptWAYsQq0w/JSVVWFmTNn4p577sE3+/Zh8+uvY/rYsWAlJWg4dw5z587F1k8/xc6//x39cnMx9Yc/RINXGf/qV7/CXXfdhb1792LgwIGYNWsWfvjDH2LJkiXYvXs3GGN46KGHPD5z7NgxrF27Fh988AHWr1+Pr7/+Gg8++KDf7C5fvhx///vf8corr+DQoUNYtGgR7rzzThQXF6tXJt5o4oZsInTbosCoPX2CBVgz015DCrHUig8zIK4LdrsrLLvrvbCJXkcayWBihHps2+ZeVZOcbHyofzlYbJsIpX1EoBUystFzFaTTydjhw4zt38/27NzJALCysjJ+rqGBsd27GfvyS8b27ePv9+9n7MsvWdvevaxz585s3bp1rksBYE888YTr/Y4dOxgA9vrrr7uOrV69miUlJbneP/XUUyw2NpYdOnTItWrpo48+YjExMayqqooxxtjcuXPZtGnTGGOMXbx4kaWkpLDt27d73Ma9997LZs6cKXmLtGrJTBih1OVMGdHUQGQiZf0TLAFCzIzKSs/z3vVFcJQ0M2bxR1CKFX3S1LLemtVpWC309CNqb+cWmeZmDEtMxHXjx2PIkCH4wYwZePX3v8d3dXWAzYbTVVW478470e+GG5A+fjzSRo1CY2MjysvLPS43dOhQ1/+ZmZkAgCFDhngcu3jxIupF0852ux05osCChYWFaG9vR0lJiU92jx07hqamJlx//fVITU11/f3973/3mLJSGxIyaqF341XS6dDUQGQRyHQP+Dr55eRwYePd4QqOkmbFqlMUYsEoLA/Ozze/T5qajt3hPuyNFrDBHGmzs/VZBZqQAAwYACQmIra1FRtfeAEf/fvfuLRHD/x51SoM+H//D6WxsZi7dCn2/u9/+NOPfoTtn36KvXv3olu3bmjxuof4+HjX/zabze+x9hADHzY2NgIA/vvf/2Lv3r2uv8OHD2vqJ0NCRk30VOpmWlZM6Esg69/o0b6WGG+ssuzeilGppaLLAvzV7D5paltvQ3UaDiZgv/99YPx4aTGjhshpbQWOHpXevqClhR/fskU/i5tIzNhaWjCqSxcsmz8fX7/7LhKSk/H+P/+Jbfv345E77sDUUaMwOCkJiTYbzp49q8rXl5eXo6qqyvV+586diImJwYABA3zSXnrppUhMTER5eTn69u3r8ZenYZ9DQkZt9IrXQlNGnhg9gtMTf9a/0aO55QXg00vbtvFXb2Fj9ikOAStOUXgPMASrV69e1hhgqGm9DXV6LZCAHTOGO04fOMAFhff3qWGlE03neIgZQcScPAnMnavvKtCEBOw6exa/feMN7D58GOXV1Vjz9dc4c+YMBuXmop/djpXFxfjm1Cns2rMHs2+7DcnJyap8dVJSEh588EHs27cPn3/+OR555BHcdtttyMrK8knbuXNnPP7441i0aBHeeustHD9+HF999RX+/Oc/46233lIlP1KQkFEbPefGacqIo9cUhJnEkpT1Tyxitm71Fbg5Ofycd4wTM2PGuCaBoAEGJ5yQD4EEbFkZj8bb2sqjDWsRp0dkAXGJmcZG/lpeDjzwAG/veq4CbWlBWn09tnz1FaY+9hj6z5iBJ5YuxR8fewxTrr0Wr7/1Fr6rr8f37rgDc5YuxSO33YaeXbqoEkumb9++uPHGG3HjjTdi4sSJGDp0KF566SW/6X/1q1/hl7/8JZYvX45BgwZh8uTJ+O9//4s+ffqEnRe/BHQFjgB0W7XEWHBvdrX2pYkSZK9I0GOVhFlXf23b5rliLSeH50dcJnY7/7v6al4Hrbhqyfs+xSuAtNrrTAWicuWdWu3R34rMnTt96/CgQbJWDSlatdTc7FoF5Prbto2xkSP17QfE+di/n69O2rePv9+9m7+XSn/4cNh7ST311FNs2LBhtNdS1CDH+Xb4cD6itJLjohXQYwrCjP4aUtY/Ye8Y8RTH1q18RLx+PTB4sLusevTQL6/hEMjKaVWH4EhGLf89f9P0I0fqE6cnIQHwtiIMHQps2KCfxU2Yzmpu5haiAQOA1FRg4ECeP8aA0lJPXx7BotSvX9TsJUVCRi2CNV67nb8vLzfPgzCS0HoKwmz+Gv5M90L9qq/3nOIQTzUKHe6aNfrkNRyCTVGUlJhPYEY7ak2vBRKwevgitrRwkSCmtBRITtZvSj8mhosRQcQkJPDjCQlczCQm8vMxXo/yaNsQUxNbkYnQdWopmIlbZNbXPJBSBBCSWT7QFIQaBAtAqAcqme5NP+0h9z537TJtuzJ9GZuVYNP0HdNLHgEHZfzWsqeWpKZzxO+bm1W8WQmcTvd3iP8X8iZMGYn/14i2tjaaWooqgjnfis36VnBctAJiB1ypEdzs2eo6WpthF/FoWXov9z4HDKB2FUnImaYXtt3QIk6Pv+kcbwdgf3Fm5OIvVk1rK/C//wFHjvDzcXFuS4yQt6NHebpos7z4gYSM3pjhQRgpiP0jdu3y7PzWrOENvKyML9lUS8yYIWJrtKyMUXKfkdquQl0pZ6YVdkoJJGBXr3avWsrPDzlOD2PM/8lA0zmCmJGazlFCoFg1zc3AhQv8uCBmAE+B5XTyvwgg4G8hExIyemOGB2GkIHbAFUZoBQXAqlXA44/zzkIQM2osiTTTLuLRsvRe7n1GYrsK1ZHZ6g7QgQTsgAHc4TY/nwelUxinR4hi29TU5P/74+K4o6xYxAio5UgbKFbNiRN80lp4f+SIe/l3c7M7T2VlqiyvNhrhtxBHGFYK2aT0xPtBuHIl72yFB6EVzeB1dVxQSD1sKip4h6LVQ1VwwBWCZMXFAX/4AzBrlqeomTUr/KkWKXO3MMUhHB83jne+ZgrQFqmI6513u/rjH4FFi6zdrgDflXLCfYjvV0gnbmOhfs5MpKdL5y09Hfj0U/dvL7ZKCFa6AH1ObGwsMjIyUFNTAwBISUlxheWX5OJF6eOCEAmH3r15v9XczMVKr148Hk5LCxAfz8WMMP105Aj/jHBcuO+mJl+xpTLt7e1oaWnBxYsXEROOFcoLxhiamppQU1ODjIwMxMbGhnwtEjJ6EYkPQmHkV1Pj+7AQOs2ePbWd5sjL4yMzQcxMn86Pi8s4SOcmC8HcDUj7awj3anW/FCsgrnerV/PAaN7CtUsXntaK7UpAEOri/kE8+PG3Ui7Uz1kFfyIHkHVPQkRaQcwYis0G1NZywSJE4I6LAzo2dMSZM0Bbmzt9bCx/L6TRIaglYwwXLlxAcnJyYNEXIhkZGZJRgpVAQkYvIvFBaJaRX14e8M9/cidPAbF/hBodtmDulrI+qSWWCHmI693ttwNdu3qKGOFh/c47XORYrV2J8R7sCHXcbpe2NAlWUH+fIwdo2Gw2ZGdno2fPnnCawc+ksZHXU4HVq/nUVVUV8NRTvluM5OTw/i47W5fsOZ1ObNmyBWPGjAlr+keK+Pj4sCwxAiRk9CLYg3DdOnc6b5RM0egx1SP+Du+Rn2DWFweo0nrk588/Qu0OO8yRYNSgdR30rncA8PzzniJGTWuc0QiOzGKhLjWt4W0FlfpcJDhAq0RsbKwqD9GwcDiAO+/k+zcJ3HmnpyjPyfEUM04nt8gkJamfH4m2Gxsbi9bWViSdO4f4jAxTtiVy9tUTf46LdXXAvfcCN90UnnOeHk5+3t/hHYju1lu5iMnPD09I1NX538VZvOrCLA64Vl4loiZ6OZqK611ZGa93UkuuI8HxWUqoV1ZyB/dAAQAj0QHaH3L6Cj2R0x8E6ruExQt2u69TsfDbq72fU6C2CwBTppjWSZyEjBlQK/y9HmH0pb4jL4872Yp5/vnQRExdHXD4MG8wU6Z4nquoAA4dcj8MDx8OHm9CjdVKcvJs5VUiaqLnVg6RuuRajNTDThwlfPRo6QjTjIUn8K0izOvr+euUKeZpe3L6g/HjuV+fVN+Vn88tboI/THm5729/+rT69+Sv7Qp+OGVl5o2SrUWkPjOha2TfcAgWyVJudFK1rqPkO9asYSwuLvxot8LGjB2bHLqidpaW8mvZ7YwlJrqvf+iQOTZy1GPTSg3QLOqsVB3Mz/dfL0Ld2NEMUZaDEFYZB6pXdrvnfYvThFsfzbpBqgQtZWXSkX2NbHve5S9sFizOU34+Y9/7Hv8dDx70/Hx5OWO5uYwlJAT+7bW4L4m269qYc8AAHkXb3z1rUB/kPr9JyISIJg8BtTpmPTp4qe+Ii+OiJlTh5LVjc0v//ryM+/fnuzoL32O3u69rlp2P9RCQKqNp+Hx/9WPnTul0Sh+MFinvsMo4mKAQtwnxdhzhChELCXNX+Yp3vzZDXRCXVWIi/63EAqS8nAsYYVd679/i0CHGrrzSGDHp1XZdA0oDxC0JmQ4sJWQYU2+vIK33HGKMsfff9/yONWv48XA6PNFnXUImOVlaxJgNC1gIxGi+D5B3HVTrwWjFB22oZexPqEtZZcRlEa7At5pQLC01X9uT+o2E/ktOXTVykCZquy4hIxaLOrU5EjIdWErIWM0iI0wXBHpIhaLSOz7rsSGcloJMTfQQkCphiEVGjQejlaY+tO4ntBQZFhDmHuVrxrbnnaecHNOKQhf+LDLbt+subknIdGAZIWNVH5n8fG6ZkfqOcEYN27ZJCxkzNnwBC3T8YnT1kRHeq+FLxZj0aFU4JlXv9JxmFOVFsozDyYve1igzigMRprfIeOfJTPmTIpCPzKBBrl3H9boPEjIdWELIqNU56dHJ6fEdHSZZDyEjNcdsJixiihejiZCRUz+0eDCayUojyktLaalnGYebFz3v0wLC3BI+MkKe/Pk0mQU/bddVhwWHam+XAg3vQ+7zm5Zfm4FAu70KS4m7dnUvNfRGWA4p5zrhRjnV+jscDr6ktLzcfSwvzx0nwm7Xb1m1XKS2nzBiGbgZCFQ/Vq3yjYkhFdcklKW/aiz7VmvJsTgvU6e6j6uxBF2vnc/NEp8pGEK/IA7AGaztab203N92NN51f+ZM85Qj4L/tChtz5ufzerVokefnzBCbSDMpZRJMY5EJ5rh18qT/80o82PVwENPqO8QjgsREt7Ov2GzcsTTbLD4QjDFzWQMUoNnUklT98J6O9DdqDqcsw7GKqf0bevt5ff658VYCuVjJmfrsWbdFRs7vpkdb9f4OcbnZ7YxlZUkvrWZM/ylQqbx7/a6ufmLHDs8wCuQjox+mEDJRtBzSA6WCR1xOBw+6Y0S0tHiW06FDphMFplkGrgDNVy0JKKm/4db1UKdDtGhj5eVu/wLBz8vsIoYx5f2VgXXfVYfLyuR9v159qdhfS3zdgwd5/JjcXN94MCYd9EhO39GqJX0xhZBRo/GIVwlJqeFdu0xV+UMWb6JO0aeMTSoIrIpuQkaqLgi/s1Rd2LUrvJFfqA6qGvg5tXz+uaeQMZtfhD/kihODrZEh1WE9/dm8y8crVpbLunzwoGkHpa4yHjNG99/ZckJm+fLlDAB79NFHXccuXLjAHnzwQda1a1fWqVMnNn36dFZdXa3ouqYQMoyF33hqa7mSl1r18ckn/qdbjHr4qyDedHvQRim6lq/4wSjVuQt1VKgfV1wReHm/P8J1UFXTwdWqFhklGGwtDrkO6+nI7C0KvaeaNmww9bSjq4zPntXd8mYpIfPFF1+w/Px8NnToUA8hs2DBApaXl8c2bdrEdu/eza6++mp2zTXXKLq2aYQMY+E1HnGHIf5budIdtt9ul24wRpkqwxRvJGS0xbDylfvwU7o6Qq2RthpLjs3iI6PHtI+BK/bCqsNGLi23wIowASP7YcsImYaGBtavXz+2ceNGNnbsWJeQqa2tZfHx8ezdd991pf3mm28YALZjxw7Z1zeVkGEsvMazc6evRSY21v2/OOqtWfxnwmiwJGS0xdDyDfbwUxqvQi3LgBoPGFFeXBYZsZ+XXm3SbEu1NRBVlrDI+MPkMXoErCBkDF9+vXDhQtxwww2YMGGCx/E9e/bA6XR6HB84cCDsdjt27NihdzbVweHgS9XEyF265nAAs2bxXVHFy/ja2vguqb168SXL48ZJ74abm6vijSggEncolrt80yo7CBuBeLn+iRPAqFHu+rpqFa/rcpb+CmUstXS0ogJIS5MfFkCtJcfivBQVSd9zuGEQ5GCmncjNtEO8GZaWy3kWRHsfoYC44Em04+2338ZXX32FL7/80udcdXU1EhISkJGR4XE8MzMT1dXVfq/Z3NyM5uZm1/v6jtgrTqcTTqdTnYx3XE/8GpRTp3hMiaoqYNAg4K9/Be6/n8c/mDSJd3jCen1vKiv5NvXCZ3/zG2D2bM80nToB/fvzhiCIv0GD+HWzsgAV710Rp04B8+cDycnuY/PnB77fDhSXsR7U1wPTpwNnzvjeg/Ab9+gBvPkmMG9e8HRr1vCHrQEYXr5ZWcBbbwETJ7qP/fGPwN13u+u6UHabNvEyE9rLRx8Bqamev8WHHwKNjfy6ZWWeZfzppzx9Sop0W/BuY4G+Nycn8H2lpLjy4uzRA4CojLOygudFLvX1/H6l8lNZyb/DO/9CvyO+z8zM8PuHYO28tpY/lKuqPPs7cb+YmMjTpaTI/lrFdVjN3zlUxPeckcFjhMXF8b5bKBvAFH0EYGw/Ifc7bYwxpnFeJHE4HBgxYgQ2btyIoUOHAgDGjRuHyy+/HC+88AJWrVqFu+++20OUAMBVV12F8ePH43e/+53kdZcuXYply5b5HF+1ahVSFDQQgiAIgiCMo6mpCbNmzUJdXR3SAog5w4TM2rVrceuttyI2NtZ1rK2tDTabDTExMdiwYQMmTJiA7777zsMq07t3bzz22GNY5B1dsAMpi0xeXh7Onj0bsCCU4nQ6sXHjRlx//fWIj48PnFjuKD6Q8q6vB44eBe65h48a8vM9rTq9egExMZ6myfx8WZYPTRBGPkJevUdgwvEAIx9FZawn3vcg/h0C3au/dAZhaPkGKhu7nVu0hg/3/ZxgaRDaiZplLMfCobAP0ayMlbavXbs8LV8ffwyMHKl/PsTHBcJoCyGVrwa/s+u6wfr5Ll0Amw04d87T8lJWxqfhWlv5ueZmU/QRgLH9RH19Pbp37x5UyBjm7FtfX88OHDjg8TdixAh25513sgMHDricfd977z3XZ44cOWJdZ99wHd0COTOKt4q3282x14jS2CESZWBqZ1+5zoJmcCr0g+lXLcl1ho3WMpa7WkjL8gnFoVhFJ1dT9RFy6/WhQ/6XY5us/jJmDWdfw1ctiRGvWmKML7+22+3s008/Zbt372aFhYWssLBQ0TVNI2TCxV+H4XC4hUxiIg+sxJg+KyS8Y4OIv0O87YLwf7DYIV4dnqk6KSnkdsgmXZ1gWPlqsZomWss42ENQj6XRSgZpKj+0TddHhFreJq2/jJGQUYy3kBEC4nXp0oWlpKSwW2+9lVVVVSm6pmFCRov4DVLXFB4KdrtbxAhoGUfGazsByb1FxBErxQHOZI7CTddJiSGLTHio2T6ivYz9PQTNtrWJFpGTzdhHKK2PJq6/jJGQMQWaCxmpPT6EKLz5+fqEczZirxOpUNuB/i8o4KHnFXRkpuykGAvNpK9zoDA5mLZ8lRDtZRzoISjX8hVow1q1+g+NRJVp67BcC4vJ6y9jJGRMgWZCJtCuq7t2uQPXicWMWYLUqYV3qO2cHHfDzcnxFDEhzNebspOS2yF/8YW5RsMS+C1fq2yAaTaLgwSG+8gE+y29p3ylrq/GoEuj4Hym7CPk9nEWqL+MWUPIGB4Qz7I0NvLXsjLPIEqHDgE/+IE7cJ1wfvt2YPRocwSpUwtxgK/ycu7xL1BZyY+JA5QJn7FygDypwGuAb7CzzEzpdGlpwOrV0kHRzBAAy0yBy4Ih97fQOvCcEVRU+Aa9vOYazyCD48bxYHf++pncXL7SUY+geenpwPr1QHGxb1vPy+PH16/n6ayMkmB7Zq6/VgvmqZOwMgzNp5bEW5tv2OC575FUmHWxQ26k4G1GDWRStbpFhjFlOwML74VXYVQq7FYuHDdgXyzJ8rXIKNGFya1HmtVhNS0cFpje8Iep+ohQ2o4Z669X3fIoY+FeRozgq6+kUDHfZJHRkro6t/WhqMituCdN4uv/BXJyfK0P3bpZf9QhRirUthhx2G0zhAZXg/T0wKNc4fdNT+ejKcHCUVLiHv3OnAkcOcKPjx8PjBmjbsj4UMnN9R3Vm2nLC2/k/hYCVhtp+kNNC0eg7SLElgKro/VvH4qFRWn91QPvrS1OneLHT51y9wMHDvD4QWax2qoim0yM6haZDrXqsRmct0UiPl7ab0Q4ZuIRjiLk+sjk5zP2n/9Ij1Z27fJczRTuqiWzjXC8R2liK53gRyW86jz6DVi+Jl9JERJ6bqLYgaksBsEw8RJgf8guX71+e7P1P6Eiav+uZ50w++C9iENDqy1ZZLRCUKtCZMrdu30tEu3t3EdC7DeSk8MjlgobO1rF+uAP8Ry93c6PVVby/+12933n5vKyuvVWvoeKeLTicHDLRJcuPIpluPPBSvw79BqZe1s4Zs0C/vAH7j/V2srTtLaab/RrdV8mKfTcRNFqhLOhrRXQ67c3o4UlFMRWJOFZV1bG32/dyv/MZLUNWzKZHE18ZMrL3Sq1c2fp0XVsrO/oRk8/A61HBnLjyHzyibtM7HZugWFMepVPuJF9lUTW1HlkLmnhMHj0G3UWGcZ09wexhEUmWnxkLHyfhrFtG2tJTuZlnJzs2U/p0EfQ8usONHP23bHD/eOKf8SdOxnLzfV9UEk95LUyMxphRpWK7Ctc33v6SGYHEtJDQE5nZZRDayCnaAM6Ur/lG+kdvo4iTVEdVmvwoeQ6VnPu9kJxHxGpAl0LOsrKQ8h4l5XG05EkZDrQPI6MWMgIVhdh/jAhga9kknqQahFpV+hsvDsnYYWMkRahEDqQgJ1UoM5aLJz8fZfeD2up+4+LY2zNGsNEQkSsWgoVnfxBdPfhUHodA/yG1CSkwY4FfYF0J5CPjPegnCwy2qOZkCkt9RUy3lFt9dr3SKozEn9nXBxjgwdLbxGgBf46R+8OZMOGgJcJGLAtWOc7eHDwzkqv0Zn3b+Ht4Ct2ANZRJEiWr8UfbLKQ+t39OeGHOeiQ/aBVS0BGyhJgmehikbFw+YSEVx1yPetKS93HvR1+NRoIkpDpQHMfmUGDPH/ExETpTlHLh4C/zmvnTvcDU08zqlR+xJaqYA+PDvx2UnI6a7n3rfXoTJzX/Hy+dYWUeBG2b9BaJIg6ZZ/yFTrlSO64vR/o3rGfpKx2YfwmhvhwRPrUoAjNyzcahL03weLI2O3uNqOx1ZaETAeqC5mOB5Nr3rC0lB/3Xoos9SNq+RAI1EiNMKPKXZodoGOV7Yzqfb9ia0egzkoPi4x3RygWCd4dodYiQU6gq0jrlMVICWDx7vHitqtSp2yYD4dZfEE0FsWaW7zMPtWqVfkGGvDs3MnYsGHS9UmY1lepHyEh04EucWQEjH4Y+PPDMKozC2SFkdERBO2kAt1vsE5HrVGrnI7ELBYOOSZjIztlrfE3uhbX08REab+2EDHUh8NoXxAdrBm6+CCZ1cKlk7XIo4yF78zPd69A9f7OK67ge3ipAAmZDjSZWqqtZS1lZf6nPYwc0UqtjDGy8W3Y4L8zDdLYZHVS3vc7eHDwhn3okDqjLCuaneU48UUy/kSllOhWoTyi2iKjgzVDt1VhZihPb3SyFnmUsc4WKhIyHWi+15KZ4kP4s1Ds3Ol7Xo+Rt5zGH6ADCckiIzVS8P4utQSI2c3O/ujIX8BlldGIBhaMqPeR0TgvsspXLYuo0RYuKXT4rX3KWMf6RUKmg6gRMuLKlZ/PLROBplm0thSoUNlD9pGRc321OjczPTSUECjQVTQid8StsN5YYtWS1mhozQhavmoNWsxokdEpb5JlrFN5kJDpICqEjFTnVVvrXgnj3XlpPf2lUmca1qolPTtrM3dyUpBFxhO5YjSEh6Lp48johUbWjKDlq0ZfYYXBiobWIr9lrIOFioRMB1EhZMzWeamUn7DiyOjdWZvR7CxFtPvIeKPkQRfCQ9H0kX31wEiLjPf3KxUiZhs0SUEWGRIyoWIqIcOY+TovFfITcmRfve/XKhaZaF+1JIVSUazwoWi6fkJvzOAj450PJW3UjIMmqTyQj0xkEzVCJgKxRBlbwewsoEYcGTMJSLVQek8KHoqWqMNKkVteZlu1FKrV1Kx1nlYtuYjRZ49tgohAKip8t66/5hr+Kt7ivqLC2HwKpKcD69cDxcVAXp7nubw8fnz9ep5Oiro6YPJkYOxYwOHwPOdw8OOTJ/N0ViI9HcjNlT6Xm+tbHnl5wMqVnsdWrvQt00hESR3o3Bno2dPdNoTyyctzt5GePXk6rXE4gDlzPI/NmeN7D1IorR96YUT5muk3FRGn67cRhJHU1QENDdKdUkUFb3xKOiWhUQPSjXrcOEMadUDS0/3fo7/OWqChAaipcQs04Z4dDregE9IZ1bnrgb+HorgORCpK6kBuLhfGUm1OEM5K21woiPNWUMBF55w5vvdgNYSBiZ7la8R3yoAsMkR0oIU1IVwLh9XIzfW1Nm3f7muVCiaIrIz3Q3HbNs/ykDPCtzJK64DR1gyrWU2VYkT5Gv2bSkBChogOvEeSwgNH/GCqqeHplGDCRq0pYhPyiRPAqFGeDwkrjmzlEukPRblYqQ6YdCqEUBcSMkR0YCVrQl2d/4dhRYXxPijR6iNCD0U3VqkD0WY1jVJIyBDRgxVGklZwqA3HcdLK0EPRjdnrgHgw4G01FQ8GBKup2QcPREBIyBDRhdlHklpNgalFtPuIRNtUohSB6sDo0dJ1QE8xoHQwYIXBAxEQEjJEdGH2kaSZp8DIR4TwVwf+8x8gMREoL+diRlwHvMWA1tYPpYMBvQYPZPXRDBIyRPRgFWuCWafAyEeE8FcH0tOBzEz+/+nT7oeytxg4dUp764fSwYAegwe9rD5RKpZIyBDRgdWsCWacAiMfEcJfHcjNBbZuBex2oLkZuPlmaTGQlqaP9UPpYEDrwYMeVp8oniIjIUNEB1azJph1Cox8RAh/dSAvj4uZQGJAz6lTpYMBLQcPety32f3rNISEDBEdWMmaYJUpMILwRo4Y0GvqVOlgQOvBg9b3bWb/Oo0hIWMkUTqfaRhWsCZYbQpML6itWAO5YkDrqVOlgwG9Bg9a33cgsfSf/wA2m/TnLN6GSMgYRRTPZxIBaG8HMjKkp8BWrwby8801BaYH1FasgRIxoKX1Q+lgQGn6cES1HlPGUmLp5ZeB+fMjtg2RkDGKKJ7PJPxQVwfcfjtw7hwXLeJRmsMBzJwJdOkCvPOOOaxHekFtxfwoEQNaWz+U+sMpSR+OqNbL6iMllu67D6isjNw2xCKcuro6BoDV1dWpet2Wlha2du1a1tLSEvpFyssZKyhgDOCv27Z5vi8vVy/DZqW2ljGHQ/JUS1lZ+GVsJRwO6d/fu574KS+lqFKH9cKibcVvGQeo98zh4OetRG0tY1dfLf1bCL/d1VczduiQqnVctfKVmz7UNqpX2w7UTux2/qewDRnZT8h9fpOQCRHVflxxxRP+TNwxq0qQzq9l0CBexmfP+n4ukh4CYnR8YFtKyDBmybYiWcZyH/pWq8dy2qXK925IHQ6ljerxm8sRS2IxI7MNWUHI0NSS0ZgxXoheBJsyKCvj7xsb3Z+JdH8JMwXDM5uDbaS0lUidKpPjTG+l1YP+CKWN6nHfcqbIcnKAV1/1/JwV25AXJGSMxqzxQvQg2HLB/HyeLifH/ZlIfQiIMcMD24yCUe22YpRQi+JlsgC0WT1YX6/vbxlKG9V61aQcsfTaa8ADD3iei4DnDQkZJajd8VG8kMCjm6Ii3/TR8BAwg7g1m2BUu60YLdTMZHmLBKZP1/e3NEMblSKQWGKMR1yOxOeNTlNdhqGaj4zXHKfHvGEoc5w6O3aanm3bPOdtt20LPDdrQX8JWZjJR0aPvMjxq9CirejU/oKWsUS9J+TjKt9Bg/TrS63oeB5GfScfmUjCe4R66hQ/fupUaCNUq4XM1xJ/oxuhjKUww/SL2pgtGJ7WwbXkWkXa29VvK2aw7Jl1VG9Fior0+S3N1kblEunPG52ElWGoumpJpF5dK2qkRgJyieTVN3IJMLpxlXG0WGR0Xs0ie6TlbTXYsEGdfCoZJardVoTrSdUju52xgwflXyNAnvyWsZqj+ijuRyQt41r2CVZecRZiPbGCRYaEjFI6KuuSuOUsO7uB9bUdZX3jSlnf/BbWty/z+zdnDmPt7epkIWII8iBrSU7mDaiszPNzVjTtykXHh5KsDsrfg14cjyIcU74Rv6X3w8hbqOXkBH8YyXygtZw961vGak5rWfnBqgI+dViPqbooE44kZEyAJnFktm1j98W+6tFe5PxFi7uLbEKJI0O+RaoRlo9MiMG1gn6PmiNpfw8ch8Odd6m4GsLxQHVIZj2UDOqopviI8vagu0UmCrGCkCEfGaV0zGsvjnsezzyzBcUJ47Et5wfYtuY0tm2D5J8w7SgOh2JJ1F61FWy54Ecf8f/T0tzHI32u1ywE8wUoL+fp7PbwV91o4e8UyP+GMf6akMDvo7ycL/EXL/MXp5NCro+N9zUBdWOKmMHXxwyIfRUjbUWOUZgtjlQgdBJWhmEGH5levXiS3bvDz4JhGGDCjprw7gYRcKQl9/fesCF8U74WI2mxpcJul7ZUxMZKW2LElqZglowgeddtNBul1ghDVi1pgdn6NFH7bykt9azDOk5ZkkVGbbxHqEKME29veQkF26kTfz1/Xrfcqo+Z4opoHViK0C+4llaxlHJz+cqqxERucRk92tNSkZsLxMYC8fGen1u9Gti6Vb5lzyyr58ySD6Po0cO6VlqjYxpJIe7vp071zI8JA4+SkJGL95RGr178eK9eQRtLRAgZMmFHH1oH19J6KWt6OpCZyf8vL3dPf9ntQEwM0NLCl3aLEZZDy53eMcsSarPkwyjWrLHutgdmGiQKiPt7YauYXbtM29+TkJFLGPPaESFkAIpGSnDUEiBa+zvl5nLrirefSmsrFzZxcUBbm7QQY0yeiDGDX4ZZ8mEkaWnWtdKadZAotENhq5iJE03b35OQUUKIUxqmEzLhOHFFuwk7UlFSJ9QSIHptIBgX5/m+spIfa20NXYjJFXOVleHlPRhWDdBGeGLWQWJeHvDXv3oeM2F/T0JGB1JT+asphEy487FWMWFbyeNeawKVBcDPKakTagoQLf2dBEuFsCpJjM3Gp5hCFWJyxZzQ+LXCbKv4qN2FjhkHiQ4HcP/9nsdM2N+TkNEBwSJjiuXX4czHWsWEbUbnOaMIVBbCFhAzZwLV1crqhNkdrsWWCrvd1yrjdHKLjPc2C3KFmFwxJw4doAV6WbXkQO0uPMw2SBTav+Aj8/HH5uzvQUJGF0w1tRTqfKyaJmytR21mdJ4zikBlIaxGqK8H3nnHfHP04SBYKux2/r683C2+hWPffitd1+QKsXDFnFrtwCyiktpd6JhtkCju7wUfmZEjTTtlSUJGB0wlZIDQ5mPVMmHrMWozq/OcGL1M8IHKQhhpFRUBV11lzjn6UElPB15/nf8viBhBfG/dysVMczNfeWVEZxyJ1gsrtDszYkY/J3F/L4QaAUy7pJ2EjA6YTsgAyudj1TJh6zVqM6vzHKD/Q8xfWQgjLSGUgBnn6MOhVy/uGyMlvpXEitECcTsYMwb44gt+3Lsd/O9/1hIzZm53ZsVsfk6AZ38v9A8CJlzSrljIzJ07F1u2bNEiLxGLKYVMKPOxapiw9Ry1mfXBbIQJXqosvFcjmG2OPlzM5D/ijdAO8vO5ZWzUKOD99z3bwapVwB13WM8yY9Z2Z1bMWk/NMmUpA8VCpq6uDhMmTEC/fv3w29/+FqcEh0HCL6YTMkbPx+o1ajPrg9kIE7xUWYhXIxhdJ7TCzJ1xXh73TRKWgk+f7iliZs2ypl+JWdudmTFzPbUAioXM2rVrcerUKTzwwAN45513kJ+fjylTpuC9996D0+nUIo+Wx1RCxizzsVqP2sz+YNbTBO+vLAQfmT17zFEnopGrrgL+9S/PYz/5iVvEWM2vxOztjohIQvKR6dGjBxYvXox9+/Zh165d6Nu3L+bMmYOcnBwsWrQIR48eVTuflsZUQsYs87EOBzB7tucx8aitooKvpgkFs4i1YOTlAS+/7HlMLObUcPoNVBaCj8y8eUBGhvF1IhpxOIDHH/c8tmCBNf1KrNLuiIgjLGffqqoqbNy4ERs3bkRsbCymTp2KAwcO4NJLL8Xzzz+vVh4tj6mEjBnmYx0O7uBYVsbN6mvWeHZ0u3Zxh9fp00O7vlnEWjAOHeKrZsQIYk4tp99AZSGsRsjM5P4ZZpujtwqBVqAFiuzrbb145RXP83/4g3VEDGCddkdEHkq31W5paWHvvfceu+GGG1h8fDwbPnw4e/nllz222V6zZg3LyMhQvme3BsjdBlwpru3jha3NA1BczHeS799f1SxYE4eDbwEPMBYXx18LChjbudPneMuAAf7LONjW9rW1PI2/PGi8/XxQyssZs9v5/QKM5eS439vt7v8LCvzfh1z8lIWrDp89G971o5naWsauvpr/TuXlnufKy1nLoEHSZSxuB971X/iLi2Ns1y797kUNdG53SvphIjSMLGO5z2/FFpns7Gzcd9996N27N7744gvs3r0bCxYsQJooguX48eORkZER9Fovv/wyhg4dirS0NKSlpaGwsBAfffSR6/zFixexcOFCdOvWDampqZgxYwZOnz6tNMuGYyqLjNGIR23CEtgTJ7hPwB/+4HZ8tNuBLl34Z7wdyuVYK8zsPCeY4MvL+X3a7e6Re04OPy6cU8M/IlBZANpHn41kgq1AE/yQvMN6i9uB2LG3oIBbKIV2cPvt1vIrMXO7IyIWxULm+eefR2VlJVasWIHLL79cMk1GRgZKS0uDXis3NxfPPPMM9uzZg927d+Paa6/FtGnTcOjQIQDAokWLsG7dOrz77rsoLi5GZWUlpoc63WAgJGREiKe2vCNFTp/OO+/8fODdd90iZerUyIoS6i3mBEFXXu4WNImJfOrHSlML0UiwFWiCH5L3Xk9CO3j7bV/H3ltv5U6ywtJs8ishiMDoZCGSTZcuXdhrr73GamtrWXx8PHv33Xdd57755hsGgO3YsUP29cwwteRwuC3FqmD2aROlbNvmaVLfto0xxlhLaSkv4+Rkbnbfts3THO9tyrcS3r+hdxls2BDadcQEqQtklleR8nLfqaGCAncd9lfGQaamWEEBP2+1Nq0TVIe1xwpTS3FBdI4P/iwiNpsNSUlJ6Nu3L2bNmoUBAwYoum5bWxveffddnD9/HoWFhdizZw+cTicmTJjgSjNw4EDY7Xbs2LEDV199teR1mpub0dzc7Hpf37Hyxel0qro8XLiWnGsmJABAPFpbgfPnnR3vQ6S+nlsuzpzhI3Zx1MVTp7j1okcPbp62wpTBqVPA/PlAcrL72Pz5QFERnD17AgCcAwYAJSWAUBcGDeL3npXFN/+zIikp/M/plC6Dxx7z/X29EdeFd97hv7cw8hfXhb/8hR/3qg9K6jARhKws4K23gIkT3cfeestdh/2VcUoK8OGHfOrJuz5nZQGffsp30BbqCuEB1WE/1Ndz6664TxCorOTnJfoEKYwsY7nfaWOMMSUXnjdvHtauXYuMjAwMHz4cAPDVV1+htrYWEydOxL59+1BWVoZNmzZh1KhRQa934MABFBYW4uLFi0hNTcWqVaswdepUrFq1CnfffbeHKAGAq666CuPHj8fvfvc7yestXboUy5Yt8zm+atUqpKSkKLlV1XA6bfjBD/jqlH/8owipqdToCIIgCCIQTU1NmDVrFurq6jz8cL1RbJHJysrCrFmz8Je//AUxMdzFpr29HY8++ig6d+6Mt99+GwsWLMBPf/pTbN26Nej1BgwYgL1796Kurg7vvfce5s6di+LiYqXZcrFkyRIsXrzY9b6+vh55eXmYOHFiwIJQitPpxMaNG3H99dcjPj4+aPq4OIbWVhuuueZ6X1+4+no+IvNWzgBXz6mpnspZGG2XlfF59L/+lUdpFd4HG8mbgcpKYMoU3zyL7s05YAA2Ll2K65cuRXxJifuzVrnHYMgoA+TnAx99JF03AJ520iS3D1F2NncUdTjc/jUOh+R1lNZhwg8B2qOrDlMZawLVYQkqK7llUOgT8vKADRv4/+K+Ii8P+Phj/31LB0aWcb3cWGJK56y6d+/OSkpKfI6XlJSwbt26McYY279/P0tPT1d6acYYY9dddx27//772aZNmxgA9t1333mct9vt7LnnnpN9PTP4yDDGWHo6nzY/csTrRKhz5H7m5C3jNyLjvltGjoxcHxnG1POP8F7K7b2c209ZkX+BCngvoxbKWajDycm8jMvKjM1nhEJ12A9S4R1yctzv7XbZ/acVfGQUr1pqbW3FkSNHfI4fOXIEbW1tAICkpCTYbDallwbArTvNzc0YPnw44uPjsWnTJte5kpISlJeXo7CwMKRrG4nflUuhbiBo9Y3ZggXme/ttoLaWv8/Pj8wooWoFJ8zLA1av9jxWWclXQVktOqzVCBYETli1lJpqUAaJqETY4d1u5+8rK90rIu12fi6C+gTFU0tz5szBvffei5///Oe48sorAQBffvklfvvb3+Kuu+4CABQXF2Pw4MFBr7VkyRJMmTIFdrsdDQ0NWLVqFTZv3owNGzYgPT0d9957LxYvXoyuXbsiLS0NDz/8MAoLC/06+poZv0JGWL4piJZx47ggmTMn8F4r/jZms9JDKz3d/0O6f38ecRbwXIYsPCDGjYuMKKGBykBu/BipuiBgJXFrRQQx2tDg+3vl5fHpvL17reF4T0QWwgDH21d19eqI6xMUC5nnn38emZmZePbZZ13B6TIzM7Fo0SL89Kc/BQBMnDgRkydPDnqtmpoa3HXXXaiqqkJ6ejqGDh2KDRs24Prrr3d9V0xMDGbMmIHm5mZMmjQJL730ktIsm4KAsWTED2dhA0HA/2jaO7S5WPiMG2ctMeOP9HS+8mrrVl9fGMFa0bkzBdgS1wW7ncfhEYfFnzkz4kZfpiOQGM3J4UKGIPTG4eDt35sI7BMUCZnW1lasWrUK8+fPxy9+8QuXI463E61dMGcF4fXXXw94PikpCStWrMCKFSuUZNOUBA2KJ0wVidWz1GhaamM2byE0bhx/0Ftlx1x/BBrFWv3e1EBcF8QmZLGgKS8HRo+OuI5Ld+rqpK0uAP8dSFQTZsLh4O2+vJy/Fxx6I7RPUOQjExcXhwULFuDixYsA4NpagAhOUCHjb6rIOzw5bcxGCAh1QRAxgk/M1q3Azp38eGIiP251fyIjqavj22GMHevbHtXa3JMg1KKiwlPE2O28PxD6BMAtZiKkT1Ds7HvVVVfh66+/1iIvEU1AIeM9VbRtm6dDq7jzNMPu1YQ5EOqCsKxaLG7z8ng92rOHxG24hOqQTxBGIAxwEhM9HXvFDsCJidwHMUL6BMU+Mg8++CB+9KMfoaKiAsOHD0cn4QndwdChQ1XLXCThV8iEMlWkhoMowbH6lIFQF6QcTnNz+R/5E4VHqA75BGEE6enAJ5/w+EZpaZ71Uhjg1Ndz38MI6RMUC5k77rgDAPDII4+4jtlsNjDGYLPZXEuwCU+E1Zc+QkZQz4D0VFGkrM4xI8KUQU2Nr4O0MNru2dNcFi5/wis9nR+vq/PNKz1gw0epQz5BGImag91AQelMMthTLGTk7GpN+CJYZBobvU4EW75Jo2nt8J4yEB5I4ikDIZ0Zyt+KwiuSkOuQTxCRQl0d389t8WJu4RHiIgGm6nMUC5nevXtrkY+IJ6CPDE0VGYPVpgysJrwijUiI3UQQSmho4JvSAnwbjg0bTNnnKHb2BYCVK1di1KhRyMnJwcmTJwEAL7zwAj744ANVMxdJBF21RBiDeLWXMGXg7a9kFgThJXYE377d18fKLMIrklDikE8QkUJuLg9ICvC9xEza5ygWMi+//DIWL16MqVOnora21uUTk5GRgRdeeEHt/FmfujqgokJayFRU0JJNM2Cl7R6sJLwiBSmH/EjbLoMg/CEEJM3PN22fo1jI/PnPf8arr76KX/ziF4iNjXUdHzFiBA4cOKBq5iyPKP5Ep+ZvAYiEDMWfMA9yY/iYBSsJr0iAYjcRBN/VXYyJ+hzFQqa0tBRXXHGFz/HExEScp3kTT0Q+DZ1+9xSADiFD8SfMgxWnDKwmvKwOxW4iCOD++z3fm6jPUSxk+vTpg70Se4esX78egwYNUiNPkYPIp6FTDXeKOl/TaLr5xajFilMGVhRekUB6uv92mptLIoaIXE6d4q9lZabtcxQLmcWLF2PhwoV45513wBjDF198gd/85jdYsmQJfvKTn2iRR2vTYX7ulM07uq9KUpF84iCSbReQfOookvvnITkZPn9ZWcCOHQbnPdKx2pSBFYUXQViJDp9GSaLRp7Gigq9WAriPjEn7HMXLr+fPn4/k5GQ88cQTaGpqwqxZs5CTk4M//elPrmB5hBd5eRj82iKk31CLOmTgIpIBBqDZ/0cuXuSbPxcW6pbL6MNqMXwoeCJBaAfFafKlc2egRw/+f1GRafscxUIGAGbPno3Zs2ejqakJjY2N6Cl0roQ0Dge6PzwTp1CNs+jOj+XZgXfece9KKuKtt4CnngKOH9c5n9GIlWL4WE14EYSVoDhNvqSn8xH11q3u1UsCJupzQoojI5CSkkIiJhiiRtCpIAu9t61G74I49HZsRe87v4/eMQ707g2Pv+HD+UeFdkMQLshXgyC0geI0SZOW5v+cSfocxULm9OnTmDNnDnJychAXF4fY2FiPP0JEiD4NBQX89cQJgDGd80wQBBGtUJwmS6J4amnevHkoLy/HL3/5S2RnZ8Nms2mRr8ggRJ+GPn34a0MDcPase4pSdeTs/CxkxKq7QxMEQSiB9tSyHIqFzNatW/H555/j8ssv1yA7EUaIPg1JSXw68tQpPhjQRMj4c2yrqwNKSoCZM4GuXQGbDfj2W57GZnPnN1qd3wiCiGxoTy3LoXhqKS8vD4zmO+QTok/DJZfwV80cfr0d2xwOLmKuvdZtTq2pAU6f5v+PHs2PT54MHDpEAf0Igog8KE6TJVEsZF544QX87Gc/Q1lZmQbZIQQEIaOZw6+UY9tnnwH79wOtrUBsLPCXv3BvdbsdKC93/02aFN3ObwRBRB4Up8myKJ5auv3229HU1IRLLrkEKSkpiI+P9zh/7tw51TIXlXT4rRQUcHHgYZFR2ydF7Ktz4gRw6638eGws0NYGTJ8O/Otfnp+prOSviYnAf/5DplaCICIDitNkWRQLGdrhWkNEfiuXLPoCQDe3RSZUn5RgDr1pab6ObS+/DDz4ILfMTJ/Oj+XkuEUMAHTrRn4xBEFEDtEcpynQc0Lc75sUxUJm7ty5WuSDADz8VgqW3wdgDbfIhBqQSU6kyvR04LvvPD/3zDPASy95bhLW0uKZJi6O1oYT1kHOCr1IfEARyrBSgEy1CPacmDIFWL4cqK/nA1gTElJAvOPHj+OJJ57AzJkzUVNTAwD46KOPcOjQIVUzF3WI/FYuqdwCgK9cujhmYmg+KVIOvYCnMDpwQHozsAcf9LzW2bP8NSfH7TNDzm+EFRA66rFjfeurw8GPT54cffvoEAQQ/Dkh+MM2NhqUweAotsgUFxdjypQpGDVqFLZs2YLf/OY36NmzJ/bt24fXX38d7733nhb5jB465mO7jR2HzqX1aEAaLitbh/h4ADG9gYnxQS/hJheIOQLEn0T+iaOYcPlK9Fz4A3T+6x9x4+mTiIuL49NH4mBPq1bxFUqCw296OiD2e/rLX4ARI9xCaNw4bnKN1NEKYX0o9DxB+EcYQIv79JUr+ZLzEyeAQYN4OontdMyCYiHzs5/9DL/+9a+xePFidBY5PV177bX4y1/+omrmopa8PNj+sRLXjNqODZiM4+gLOAEcC+Vi8QD64gj6Yv25KcCvAOAV/LFLLhb3eR+orXXHiDl8GJg1i4uYuDigb1+gtNTzcrfdxi035PxGWIVgHXW0rr6j6TZCwHvhh+AzWVDAN4vct8/I3AVFsZA5cOAAVq1a5XO8Z8+eOCtMPxDh0RGQaQ2q8BW+h3bEANk5wIsvur3qFdDWBnz9dgk+/+thnEIvfImr8Fr6Yiza9DBsjQ3c12XsWB4AT5gDfeEF4JFHgOZmPpX0wgtcxLS2ArffDmzZEtnOb0RkEaijjsZAZ7TTM+GNv4jGvXpFnpDJyMhAVVUV+ghx9Dv4+uuv0ct7d0xCOSJzd0pBAUavfLZj5Pg58NMvQ+t0HQ6MXz8Ri1GOenRGFqrxTVkKdm89jSuHMU/z+rp1wPnzwB13uH1nhO/cto2LmLIy95QSdXLWIdpH4BR63g1NtxHe+ItovGmTMflRgGJn3zvuuAM//elPUV1dDZvNhvb2dmzbtg2PP/447rrrLi3yGD1oEZBJ6JjKywG7HWk5nXEr3gcAvHn7R9wfRvx9l14K9O/PR2Peo9WrruKWmIICmlKyGoEcXk+dig6HV38ddTQ6rNNOz4SYQBGNp041OndBUWyR+e1vf4uFCxciLy8PbW1tuPTSS9HW1oZZs2bhiSee0CKP0YPaAZmkhBGAucOXYNWZ2fhn0y1oL78ApKYC19wK/Da144PpwKWbgb5O0TEgORl4+OE89KEpJeshNQLPyuLnpk6N/BG4d0ct9pERWySiCZpuIwDp54S4bgirliorgd69DctmIBQLmYSEBLz66qt48sknceDAATQ2NuKKK65Av379tMhfdKF2QCY/wui6fz+I3DEOVCAPr+ABoBHAP7w/nNjx50lTE/DKKzRKsxxSDq9vvcXPiacQI3EEHqyjjubVdzTdRgQbQE+axN+npkp92hQoFjICeXl5yKPKrj5qBmSSEkYOB2LnzcE6dMY63AQGG9ClC3D3PQEF0r59wJo10WmFjxi8H94TJwKrVwP5+cCGDZH78KLQ8/6hnZ6JYAPojz4C9u7lUeBNSshChrAIYmEkMq9fXlCAy1dOcZvX1/4lYOf14YdcyJw+rVvOCS2QGoH/9a+R/dCK5tDzgaDpNkIg0AA6J4cLGRMTUmRfwoKE6UicmclfSchYHKkR+P33R76pLT3dv0UzNzf6RAzt9ExEECRkogXBvO7tyCeY14OsRBKETE0NbbFkWbxH4B9/zI8Ly+kjXcwQbsLsDwjCTNDUUrQQpnldcDFoaeHBgLt00Ta7hMpIjcCzsnjUzvx84JtvotfhNRqh6TYiggjLIjNkyBA4aBRnHcIwrycluU/T9BJ4vBV/ZveKCvPFY/E3Age4mKERePRB021EhBCWRaasrAxOp1OtvBAmJzOTP59raoCBA43OjYFYMbx7oBF4r140AicIwrKQjwwhG2F6KeotMsG2vT9xgp9vaDAyl77QCJwgiAgkLCHz/e9/H8nJyWrlhTA5tHKpAy3Cu1ttqoogCMIkhCVkioqKkJ2drVZeCJNDQkaEeHWHEN7dO3KsXALtg+RwRMc+SARBECFCU0uEbEjIeCEElxMTSnh3q05VEQRBmAASMoRsSMh4odZuyrQTMUEQRMiQkCFko6mQCeQjUllpvmmVQNvehxJcTs2pKoIgiCiChAwhG82ETCAfEQCYMsVcPiJahXdXa6qKIAgiighLyDQ2NqK+vt7jj4hAOqwl4m0KXKixosafj8ipU/y1rMxcPiJahXdXa6qKIAgiilAsZEpLS3HDDTegU6dOSE9PR5cuXdClSxdkZGSgC8WtjzxE1pKeLdzC0NQENDZCvRU1/nxEpk7l5/Pz9fMRkbMMWgguV1zsay0RwrsrDYan9lQVQRBElKA4su+dd94Jxhj+9re/ITMzEzabTYt8EWZBZC1JvWEsUpKPoulCDE5/XYnUeeP4g1ZIF05ANcGaITzMR40ChBhFRUX6TK8ojdjr736VCi6pqSrv8qB9kAiCICRRLGT27duHPXv2YMCAAVrkhzAbgrWk44GaGVeBUthx+o5HcUmlyitqBB+RUaM8j/fqFf615eA9xSUICrG1REinZhRcYaoKkJ6qEgQU7YNEEAThg2Ihc+WVV8LhcJCQiSZED9TME6dQCjvuqlyOtISngJT+wLQEn4/07Qv84x9Agu8p/0j5iADcVyY/P9Tcy8dLtGHcOC6s5szRdhk07UQcGdTVSf+GAF95l5FBvyFBaIBiIfPaa69hwYIFOHXqFC677DLEx8d7nB86dKhqmSNMRIe1ZOio/diJQhxHX6AFwEHp5F9/DTz0EDBmjMzre/uIrFwJzJ/Pz02dCmzYoM/0ktQUF6D9Mmg1p6oI/Qk0LQnwlXeCYCUxQxCqoljInDlzBsePH8fdd9/tOmaz2cAYg81mQ1tbm6oZJExCh7XkRVTgdrwDJ+KBrGzg978HevTwSPr448DBg16rmwLhz0ekqAjYt4+vWtLTR0RqiouWQROB8DctKV55l52t/rQkQRDKVy3dc889uOKKK7Bjxw6cOHECpaWlHq9EBCKyliQW5OLabb/GpIJjmFT9FiY9dTUmXerApElw/fXtyz925ozM6/tbziz4xuTn6+sjQsugCaWYaeUdQUQZii0yJ0+exH/+8x/0FZ5WRGQTwooawUAjW8gE8hEBgI8+0s+/QGqKS/CREY+0CcIbM6y8I4goRLFF5tprr8W+ffu0yAthRkII/qZYyABcpPgbrebk6CNitIrYS0QPUtGZAf1W3hFEFKLYInPTTTdh0aJFOHDgAIYMGeLj7HvzzTerljnCBISwoiYkIWMGaBk0ES5Gr7wjiChEsZBZsGABAODpp5/2OUfOvhGKwhU1lhUytAyaCAezrLwjiChD8dRSe3u73z8SMQQAdO/OXy0nZAAuUjp3lp4+ys3lIscsm1cS5sHftGRRET8vrLyjaUmCUB3a/ZpQHctaZIDAO3GrtbcUEXmYbeUdQUQRsqaWXnzxRdx///1ISkrCiy++GDDtI488okrGCOsiCJmzZwHGANW24woUObWiQp1pH6O2KSCsjZlW3hFElCFLyDz//POYPXs2kpKS8Pzzz/tNZ7PZSMgQLiHT2sq1R0aGChdVuqFjqBi1TQFhfQL5kuXkAF4LIwiCUAdZU0ulpaXo1q2b639/f0oD4i1fvhxXXnklOnfujJ49e+KWW25BSUmJR5qLFy9i4cKF6NatG1JTUzFjxgycPn1a0fcQ+pKUBKSm8v9Vm17ytpQI0z5iS0lNDU8nUFfn3yehosL/9JB4abkQD8Q7jg5BEARhChT7yDz99NNoamryOX7hwgXJlUyBKC4uxsKFC7Fz505s3LgRTqcTEydOxPnz511pFi1ahHXr1uHdd99FcXExKisrMX36dKXZJnRGdT8Zf5FTvR0sBUtJuL4ueXnAyy97HhNvUxBICBEEQRC6oVjILFu2DI2NjT7Hm5qasGzZMkXXWr9+PebNm4fBgwdj2LBhePPNN1FeXo49e/YAAOrq6vD666/jueeew7XXXovhw4fjjTfewPbt27Fz506lWSd0RBOHXyWWklAsOGIOHQK8YyIJ2xSQ0y9BEIRpUCxkhM0hvdm3bx+6du0aVmbqOh4KwnX27NkDp9OJCRMmuNIMHDgQdrsdO3bsCOu7CG3RbOWSVORUqQ0dQ7HgCNNQDgeP+9HczN9nZvJ0J04Ao0fzv2BCiCAIgtAF2QHxunTpApvNBpvNhv79+3uImba2NjQ2NrqC5YVCe3s7HnvsMYwaNQqXXXYZAKC6uhoJCQnI8PIWzczMRHV1teR1mpub0Sw8gADU19cDAJxOJ5xOZ8j580a4lprXjCS6dYsFEIPq6jY4ne0hXUOyjE+d4kHGhD1sAP6+qMg3DHxWFrBpExclZWWAIIgHDeLps7IApxOorwemT+eq629/A+65h//fvz/Q1gZUVgIJCXwJreCf1b8/D3CWmcmvYUGoDmsPlbG2UPlqj5FlLPc7bYwxJifhW2+9BcYY7rnnHrzwwgtIF3nnJyQkID8/H4WFhaHlFsADDzyAjz76CFu3bkVuxyh51apVuPvuuz2ECQBcddVVGD9+PH73u9/5XGfp0qWSU1yrVq1CSkpKyPkjlPHmm5di7dp+uPnmY7jnnkNGZ4cgCIKwGE1NTZg1axbq6uqQlpbmN51si8zcuXMBAH369ME111zjs8dSODz00EP48MMPsWXLFpeIAYCsrCy0tLSgtrbWwypz+vRpZGVlSV5ryZIlWLx4set9fX098vLyMHHixIAFoRSn04mNGzfi+uuvV7UsLE99PdDYiG9G5mLtWiA1tQBTp/bm5yor+XImmb+DRxmfOQNMmcItK/n5bgvMqVNui0t+Po/XkZPjvoj4vID481Lp7HbgmWeAn/+cv8/OBqqq3Gnffx+49lrlZWMyqA5rD5WxtlD5ao+RZSzMqARD8V5LY8eORVtbG9577z188803AIBLL70U06ZNQ1ycsssxxvDwww/j/fffx+bNm9GnTx+P88OHD0d8fDw2bdqEGTNmAABKSkpQXl7u1/qTmJiIxMREn+Px8fGa/AhaXdeS1NUBN94I1NQga+EXALrh229jEB8fw/1Orr02pFgv8fHxiBeCiWVne+5Zk5/P348bx89nZLjjdTgcwHXXee59M2cO8M03/LjYQVh8nZIS4NZb+XG7nQfEuXDBnaEHH4yoZdhUh7WHylhbqHy1x4gylvt9ioXMoUOHcPPNN6O6uhoDBgwAAPzud79Djx49sG7dOpd/ixwWLlyIVatW4YMPPkDnzp1dfi/p6elITk5Geno67r33XixevBhdu3ZFWloaHn74YRQWFuLqq69WmnVCa0QrhXr8/icAXufOvmpExVW6oaPU3jfiXawFB+DiYvf1BEfiUaPc125t5ZYksRDyjvpLEARBGIZiITN//nwMHjwYu3fvRpcuXQAA3333HebNm4f7778f27dvl32tlzvidIwbN87j+BtvvIF58+YB4FGFY2JiMGPGDDQ3N2PSpEl46aWXlGab0ANRVNweJw4AAKpOtmBP4UPAqQyg183Ai/8HnM4CZMQ0bG0Fjh1Lx1dfAdzYl87/JD+b63m8MR1I+T7Q6zL+nTVZQA0A5AEvbgN++EMgpQtwXHS96mrgh78D8D33dSoBZF0OvPg6ehVkISuQECIIgiB0R7GQ2bt3r4eIAfiKpt/85je48sorFV1Ljp9xUlISVqxYgRUrVijNKmEEHVaPHqPuBBxA1bcJGIEP+LlTAG5UcrF4AONCzEhnAG/yf32+MwvABzw/4ySOe1PNrxEbC5SU5OESQczQJoAEQRCGo1jI9O/fH6dPn8bgwYM9jtfU1KBv376qZYywMHl5yF+9HD8Y/S/sQIcvU8+eQIKv71JgGC5cuIjk5CQAau08KUFbK1Bzhr/GxgE9e/BXpxM4ewZoa8NpZMLZloCSEuCSqRJTWQRBEIQhKBYyy5cvxyOPPIKlS5e6/FR27tyJp59+Gr/73e88vIzVXCVEWAiHA7a75uBfEO29lap8nyKnsxVFRR9j6tSp2jqZ1Z0HJk8XbUgpNIt4wNEGjBuH0WfWYFvDMFy82HGKppMIgiBMgWIhc+ON3E5/2223uYLiCVNEN910k+u9zWZDW1ubWvkkrILYsdcqDrIyHImT7swCiuEWMgRBEIQpUCxkPvvsMy3yQUQCoawUMgvp6f6niXJzkdSxmzcJGYIgCHMRUhwZgpCkc2fuCwN4Wl7EYsaiDrJJSfyVhAxBEIS5UCxkAKC2thavv/66KyDe4MGDcc8993hsW0BEIUpjvVgIEjIEQRDmRPHu17t378Yll1yC559/HufOncO5c+fw3HPP4ZJLLsFXX32lRR4JK5Ge7n/aKDfXkiIGICFDEARhVhRbZBYtWoSbb74Zr776qmtLgtbWVsyfPx+PPfYYtmzZonomCcJoSMgQBEGYE8VCZvfu3R4iBgDi4uLwk5/8BCNGjFA1cwRhFkjIEARBmBPFU0tpaWkoLy/3Oe5wONDZgk6cBCEHEjIEQRDmRLGQuf3223HvvffinXfegcPhgMPhwNtvv4358+dj5syZWuSRIAyHhAxBEIQ5UTy19Ic//AE2mw133XUXWltbAfCtth944AE888wzqmeQIMwACRmCIAhzoljIJCQk4E9/+hOWL1+O48ePAwAuueQSpKSkqJ45gjALJGQIgiDMSUhxZAAgJSUFQ4YMUTMvBGFaSMgQBEGYE8U+MgQRjZCQIQiCMCckZAhCBiRkCIIgzAkJGYKQAQkZgiAIc0JChiBkQEKGIAjCnJCQIQgZkJAhCIIwJyRkCEIGJGQIgiDMCQkZgpABCRmCIAhzQkKGIGQgCJnmZmPzQRAEQXhCQoYgZEAWGYIgCHNCQoYgZCAImZYWoL3d2LwQBEEQbkjIEIQMBCED0PQSQRCEmSAhQxAyEAsZml4iCIIwDyRkCEIGcXFAbCz/n4QMQRCEeSAhQxAyIYdfgiAI80FChiBkQkKGUEx9vf9zFRVAXZ1+eSGICIWEDEHIhIQMoYi6OmD6dP7/qVOe5xwOYOxYYPJkEjMEESYkZAhCJiRkCEU0NABnzvD/p07l4gXgr+PGASdOADU1PB1BECFDQoYgZEJChlBEbi5QVMT/Lyvj4mX7dreIKSgANm/m6QiCCBkSMgQhExIyhGJ69eKv+flcvIwa5Sli8vKMzB1BRAQkZAhCJiRkiJD56189369cSSKGIFSChAxByISEDBEy99/v+X7OHLfPDEEQYUFChiBkQkKGUIywWqmsjE8nbdvGX0+c4L4yJGYIImxIyBCETEjIEIqoqOCrlQDuI7N5M3DNNfxVLGYqKozLI0FEACRkCEImJGQIRXTuDPTowf8vKnL7xOTlucVMz548HUEQIUNChiBkkpjIX0nIELJITwfWrOH/C6uXBPLygOJiYP16no4giJAhIUMQMjHMIlNX53/6gcLcm5u0NP/ncnNJxBCECpCQIQiZGCJk6up4GPuxY30dQynMPUEQBAkZgpCLIUKmoYGHsfde5UJh7gmCIACQkCEI2RgiZHJzfVe5UJh7giAIF3FGZ4AgrIJhPjLCKhdBvIwaxY9TmHuCIAiyyBCEXAxdfp2Xx8Pai6Ew9wRBECRkCEIuhgoZh4OHtRdDYe4JgiBIyBCEXAwTMmLHXgpzTxAE4QEJGYKQiSFCpqLC17GXwtwTBEG4IGdfgpCJIUKmc2cexh7wdOwVOwBTmHuCIKIYEjIEIRNDhEx6Og9j39Dgu8RaCHPfuTNFiCUIImohIUMQMjHMRyY93b9QofgxBEFEOeQjQxAyod2vCYIgzAcJGYIIRsemjZJChjZtJAiCMBQSMgQRCNGmjUnfVQEQCRnatJEgCMJwSMgQRCBEmzYmzfkBgA4hQ5s2EgRBmAJy9iWIQAibNo4bh6QTJQCAlhbg50PXwVZ7D5DRBZh6J/BymusjNhswbRpw5ZUG5ZkgCCKKICFDEMHoiNmSOmYS4sta4EQCltc+yM/VAviL70fefRcoKdEzkwRBENEJCRmCkENeHpL/+RpWj5qJz/F9fuwHtwE5OR7JmpqAV18FSksBxrh1hiAIgtAOEjIEIYeOTRtn4ARmYA0/tufPwB83e+xA3dLChYzTCXz7LdC9uzHZJQiCiBbI2ZcggqFg08aEBKBbN/5/VZUhuSUIgogqSMgQRCBC2LQxO5u/VlYakF+CIIgog4QMQQRC2LRREDHemzYWFPhs2igIGdUtMh2B+SShwHwEQUQphgqZLVu24KabbkJOTg5sNhvWrl3rcZ4xhieffBLZ2dlITk7GhAkTcPToUWMyS0QnwqaNxcUevjAA3Js2rl/vsReS4P+rqpARBeYTT2UBoMB8BEFENYYKmfPnz2PYsGFYsWKF5Plnn30WL774Il555RXs2rULnTp1wqRJk3CRNrsh9CQ93f/mjLm5Phs6amKREQXm8/DLocB8BEFEOYauWpoyZQqmTJkieY4xhhdeeAFPPPEEpk2bBgD4+9//jszMTKxduxZ33HGHnlklCNloImREgflcYmblSmDOHE//HdoNmyCIKMO0y69LS0tRXV2NCRMmuI6lp6dj5MiR2LFjh18h09zcjObmZtf7+vp6AIDT6YTT6VQtf8K11Lwm4YlVy7hHDxuAOFRWtsPpbFPvwllZwKZNwNSpQFkZILSNQYOAoiJ+XkFZWbV8rQSVsbZQ+WqPkWUs9ztNK2Sqq6sBAJmZmR7HMzMzXeekWL58OZYtW+Zz/OOPP0ZKSoq6mQSwceNG1a9JeGK1Mi4r6wrg+zh+/AKKij5R/wuWL/c9tm8f/wsBq5WvFaEy1hYqX+0xooybmppkpTOtkAmVJUuWYPHixa739fX1yMvLw8SJE5GWlhbgk8pwOp3YuHEjrr/+esTHx6t2XcKNVct4wADg5z8H6utTMGXKVHWj+5465bbICOTnc4tMr16KLmXV8rUSVMbaQuWrPUaWsTCjEgzTCpmsrCwAwOnTp5EtOB10vL/88sv9fi4xMRGJiYk+x+Pj4zX5EbS6LuHGamVst/PXCxdsaGqKR0aGShd2OIDrrnP7xAg+Mt98w4+Ll4crwGrla0WojLWFyld7jChjud9n2jgyffr0QVZWFjZt2uQ6Vl9fj127dqGwsNDAnBFEYFJS3AuZVHP4DSEwH0EQRDRgqEWmsbERx44dc70vLS3F3r170bVrV9jtdjz22GP49a9/jX79+qFPnz745S9/iZycHNxyyy3GZZogAlFXBzQ0IDs7F3V1XMgMGtRxrqKCB87zWq4tCyEwHyAdmG/cOJ/AfARBENGAoUJm9+7dGD9+vOu94Nsyd+5cvPnmm/jJT36C8+fP4/7770dtbS1Gjx6N9evXIykpyagsE4R/hKB1NTXIzjqEI0hyW2SEeC89e/oE0JOFEJivocF3ibUQmC9UkUQQBGFhDBUy48aNA2PM73mbzYann34aTz/9tI65IogQEQWty67ZAGAaFzLioHVCulAER3q6/89R/BiCIKIU0zr7EoTlEAWtyznxPwDA+3+vR91vPgBq5wEZXYAb5wD/Zw6rSVtbDI4dG4hdu2IQG2t0bpSRkwP88IdAjGm9/AiC0AsSMgShJh0+K/lXvAZ8C2w/kIbteIifqwXwopGZ8yYWwACjMxEyubnATTcZnQuCIIyGhAxBqE1eHuasnoozE5fiHLryY//vB+69C0xCe3sbyspOIj+/N2JirGOS2bED2L0b+PRTEjIEQZCQIQj1cTiQtmAWluKE+9hXfwo5zotWOJ3tKCo6gKlT8xAfbx0h8/bbwMyZ3L+ZIAiCZpgJQk3Ejr0FBcC2bZ5xXoRdq4nQqKvD2H6VAIC9e4HaWtG5igq+cowgiKiChAxBqAUFrdOWjuXt2bd9H/36OMEY14kAuEAcO5YvfycxQxBRBQkZglALIWidIGK8g9YVFFDQunAQLW8fe3YNgI7pJbEVrKaGpyMIImogHxmCUAsKWqctouXtY058iNdwOz7693mMXbkUqB4IZI0DfvkMsK8HENpG4KrT2mrD7t2ZAGyIi+PVYOhQo3NFEJEFCRmCUBMKWqctHdatsaNnAuXAwROdcCNe5+eqAdxtaO4kiANwtceRgweBwYONyQ1BRCIkZAiCsBZ5ebCvfhY/G7Ucm3AdPzZoEJBqvik7xtpRV1eH9PR0/O9/MaivB/bvJyFDEGpCQoYgCGvhcABz5mC5eHl7cwGwYbOplrcDgNPZhqKiLZg6dSruuisGb7+t4o7oBEEAIGdfgiCshIWXt+fk8NfKSmPzQRCRBgkZgiCsgcWXtwtChiwyBKEuNLVEEIQ1EJa3A9LL28eNM/XydmGHCrLIEIS6kJAhCMIaWHx5O1lkCEIbSMgQBGEdLLy8nSwyBKEN5CNDEAShA4JFpqEBaGw0Ni8EEUmQkCEIgtCBzp2BTp34/yFPL9XV+Xdmpk0ziSiFhAxBEIROhLUEu2PTTIwd67vMnDbNJKIYEjIEQRA6IfjJhGSREW2a6REzhzbNJKIcEjIEQRA6EZZFRtg0UxwzZ/t239g6Jnd6Jgi1ISFDEAShNfX1QEWF9BJsJb4tQswcQcyMGuUpYky2RQNB6AEJGYIgCK2ZPh0YOxbZybUARBaZUHxb8vKAlSs9j61cSSKGiFpIyBAEQWjNmTPAiRPIeXUpgA6LTKi+LR2bZnowZ46p95kiCC0hIUMQBKE1RUVAQQGya/YBACqPN4Xm22LhTTMJQisosi9BEITW9OoFbN6MnGvuBiqAk+U2PIJHgbR0YMx04Pcy9odqbAD+/SlQL/rc252BMXuBs2uAE3XAZZ8CM6YDqcGvl5QEPPAA0KdP+LdHEEZCQoYgCEIP8vLQ641fI+H6ZlxEMv6MR4B6AG/KvUBnAHP5vx6f8zr+hvwsffst8Prr8tMThBkhIUMQBKEHDgdSfzgba9Ef23ENP5bRBbhrDrewyOHiRaClBUhL8z1XXw8kJHBTSxAOHwbWrAFOnlSQf4IwKSRkCIIgtObUKeC664ATJzClAJiy8pfcQffECeDD5xUsnU7q+JNCQtz44ZNPuJCprpb9EYIwLeTsSxAEoTVTp3o69l5zjW9wO397KGlAWBGGCcJkkJAhCILQmh49fIPWiYPb9ezJd5XUiaws/nruHNDcrNvXEoQm0NQSQRCE1qxZA1y44LvEOi8PKC7mIiZdpp+MCnTpAsTHA04nD2FDsfQIK0MWGYIgCK1JS/MfJyY3V1cRAwAxMUBmJv+f/GQIq0NChiAIIgoRppfIT4awOiRkCIIgohDB4ZcsMoTVISFDEAQRhQgWGRIyhNUhIUMQBBGFkJAhIgUSMgRBEFEI+cgQkQIJGYIgiCiELDJEpEBChiAIIgohZ18iUiAhQxAEEU3U1QEVFR5TS4x1nKuo4OcJwkKQkCEIgogW6uqAyZOBsWOR6eR7OzU3d2gXhwMYO5afJzFDWAgSMgRBENFCQwPfk+DECaRMGYu0zu0AgOqvq/jGlSdO8PMNDcbmkyAUQEKGIAgiWsjN9dh1O/tiKQCgevaPPHfn9redAkGYENo0kiAIIpoQdt0eNw5ZJxwowSWYVvUyEmP+DNRlAMNjFV0uOxvYsMG9Coog9IaEDEEQRLSRlwesXIlrRm1GMcahHulAO4BvlV/qzBngv/8F7r1X9VwShCxIyBAEQUQbDgcwZw5+gxOYhzfRggQgNw944w1FppVnnwVWrgSOHNEwrwQRBBIyBEEQ0YTD4XLstRUUoP/KN4E5c4AT64EfjuLTTnl5si519dVcyJSUaJlhgggMOfsSBEFECxUV7tVJgmPvNdd4OABj3DieTgYDB/JXssgQRkJChiAIIlro3Bno2dMtYgTLi+AAXFDAz3fuLOtygpA5cYLHoyEII6CpJYIgiGghPR1Yv57HifFeYp2XBxQXcxGTni7rctnZPHlDA3D8OHDppRrkmSCCQBYZgiCIaCI93X+cmNxc2SIGdXWwnaqQnl6irQ4IHSEhQxAEQShDtNXBwLzzAERChrY6IHSGhAxBEAShDNFWBwM+ewVAh5ARrYiirQ4IvSAhQxAEQShDtNXBwO+2AwBK9jT4roiirQ4IHSBnX4IgCEI5HSudBl5zH1ABHDwcg0fxCJCWDoydAfxB3sqncGhvj0FZ2WX45JMYxNCwXBPklvHtt/OV/EZAQoYgCIIIjbw89F35FJLHN6EJnfAiHgXqAbyhVwZiAVyi15dFKfLKePBgEjIEQRCE1XA4kHjvnfgP+uAzjOfHMroAd94JpKVp/vVtbW04fvw4LrnkEsTGKtvskpCH3DK+4godM+UFCRmCIAhCOSLH3gkFwISVT3dsdXACKPqjoq0OQsXpbEdR0TeYOrUP4uNJyGiBFcqYZhUJgiAIZai81QFBhANZZAiCIAhlCFsdANJbHYwbp2irA4IIB0tYZFasWIH8/HwkJSVh5MiR+OKLL4zOEkEQRPQibHVQXOw7fSRsdbB+vfwowQQRBqYXMu+88w4WL16Mp556Cl999RWGDRuGSZMmoaamxuisEQRBRC9qbXVAEGFi+qml5557Dvfddx/uvvtuAMArr7yC//73v/jb3/6Gn/3sZ7Kvc77lPGJbfB2VYmNikRSX5JHOHzG2GCTHJ7veX2y7iPMt5xHP4oOmbXI2gTEmeV2bzYaU+JSQ0l5wXkA7a/eb504JnUJKe7H1Itra21RJmxKfApvNBgBobm1Ga3ur7LSByjg5PhkxNq7FW9pa4Gxz+r2ukrRJcUmIjYlVnNbZ5kRLW4vftIlxiYiLiVOctrW9Fc2t/rcWTohNQHxsvOK0be1tAcs3PjYeCbEJ7rStF/1eV5y2nbXjgvOCKmnjYuKQGJcIAGCMocnZpEpaJe0+nD7ifMt5v2VMfYQb6iM4ZuwjAtVhPfoIOdiYv9ZgAlpaWpCSkoL33nsPt9xyi+v43LlzUVtbiw8++MDnM83NzWgW7SdfX1+PvLw84GcAknySY8olU/DB7e7rZPw+w28HOMY+Bp/c+QkAwOl0Ivu5bNS31UumHZ49HDvu3uF6329FP5ysOymZdlD3Qdh3/z7X+2F/HYZvzn4jmbZ3em8cXXjU9b7wjULsqdojmbZ7cndULqp0vZ/wjwnYUr5FMm1KfApqf1zrej/tnWn46PhHkmkBoOXn7gZ2x5o7sObIGr9pv3v8O1endu+6e7HywEq/aU89ego9OvUAADxU9BD+uvevftP+78H/IT8jHwDws00/w3O7nvOb9uv7vsbgHoMBAE9veRq/3vprv2m3z9uOETkjAAB/3PlHLPl0id+0G2dvxNjeYwEAL+9+GY9+/KjftGtvW4upfacCAP6+/++Y/+F8v2lX3boK/2/Q/wMAvPfNe5j1/iy/aV+78TXcNfQuAEDRsSLc8q9b/Kb908Q/4YERDwAANh3fhCnvTPGbdvm1y/Gjq38EANhduRvXvOk/SMQTo5/Ak2OeBAAcOnMIV7zqfy3m4pGL8cx1zwAAymrL0P+l/n7TLvjeArw4+UUAwJnzZ9DrT738pp0zZA5ev+l1AFxAdPlDF79ppw+cjrenv+16n/Bb/51mqH0EAOQ8n4OzF85KpqU+wg31ERyz9RHFJ4tx/T+v95tW6z6ivr4e3bt3R11dHdICLOc3tUXm7NmzaGtrQ2ZmpsfxzMxMHPHYatXN8uXLsWzZMtnfUXOmBkVFRa73bW3+Rw3ffvutR9pA1NXWeaRtavI/OmxsbPRI29jY6DdtU1OTR9q6Wv+bsrW0tHik/fbbb/2mbWtr80hbcybw1J04bXVVdcC0GzZsQFIsV5EVQVYxfLLpE6THcZO0o8IRMO1nn32GzEReN06cOhEw7edbPsfJZP6QOFp1NGDabdu3oSaF3/+RGul6JrBz506cP8RH6IfOHAqYdveXu4H/8f/3fbsvYNqvv/4aKaV8VP117dcB0+7btw9FFfz32F23O2DaQ4cOoaiGpz3QcCBg2iNHjqDoHE97tClwmR09ehRFjTxt+YXygGlPnDiBomae9nTz6YBpT5afdNW1utbAGxBWVFS40l5s8z8yBHidlduWw+kjWlr8j6ipj3BDfQSH+giO0EcEahNiTG2RqaysRK9evbB9+3YUFha6jv/kJz9BcXExdu3a5fMZfxaZk5UnJRVdqGZjp9OJdevX4dprr0V8PE0taWE2brzQiI8/+dhvGZPZmBOq2fhi80UUfVzkt3xpakl5Wu92X3u+Fp9++qlkGVMf4Yb6CI7Z+oi29jY0XGjwW4e17iMiwiLTvXt3xMbG4vRpz1Hb6dOnkZWVJfmZxMREJCYm+hzP6JSBtE7BI01mxGfIzl9SbBIyOmVINiBv0uPlO74pSSvnu62aNhWpssvYDPmNj49HClKCJ1SaFvFITkwOnlBhWkB+HY5HPJISJeZm/ZCY4NsG1UibkCB/3lxJWiXtXlHaThmyy5j6COVpqY/oSKtRHxGPeC7kZZSxFn2E3LI19aqlhIQEDB8+HJs2bXIda29vx6ZNmzwsNARBEARBRCemtsgAwOLFizF37lyMGDECV111FV544QWcP3/etYqJIAiCIIjoxfRC5vbbb8eZM2fw5JNPorq6GpdffjnWr1/v4wBMEARBEET0YXohAwAPPfQQHnroIaOzQRAEQRCEybCEkAkHwbO/vl463kuoOJ1ONDU1ob6+XpGzFyEfKmNtofLVHipjbaHy1R4jy1h4bgdbXB3xQqahoQEAeFA8giAIgiAsRUNDA9IDbHlh6jgyatDe3o7Kykp07tzZFXtADYT4NA6HI+D6diJ0qIy1hcpXe6iMtYXKV3uMLGPGGBoaGpCTk4OYGP+LrCPeIhMTE4NcfxubqUBaWho1II2hMtYWKl/toTLWFipf7TGqjANZYgRMHUeGIAiCIAgiECRkCIIgCIKwLCRkQiQxMRFPPfWU5HYIhDpQGWsLla/2UBlrC5Wv9lihjCPe2ZcgCIIgiMiFLDIEQRAEQVgWEjIEQRAEQVgWEjIEQRAEQVgWEjIEQRAEQVgWEjIhsmLFCuTn5yMpKQkjR47EF198YXSWLMnSpUths9k8/gYOHOg6f/HiRSxcuBDdunVDamoqZsyYgdOnTxuYY/OzZcsW3HTTTcjJyYHNZsPatWs9zjPG8OSTTyI7OxvJycmYMGECjh496pHm3LlzmD17NtLS0pCRkYF7770XjY2NOt6FeQlWvvPmzfOp05MnT/ZIQ+Xrn+XLl+PKK69E586d0bNnT9xyyy0oKSnxSCOnXygvL8cNN9yAlJQU9OzZEz/+8Y/R2tqq562YEjnlO27cOJ86vGDBAo80ZipfEjIh8M4772Dx4sV46qmn8NVXX2HYsGGYNGkSampqjM6aJRk8eDCqqqpcf1u3bnWdW7RoEdatW4d3330XxcXFqKysxPTp0w3Mrfk5f/48hg0bhhUrVkief/bZZ/Hiiy/ilVdewa5du9CpUydMmjQJFy9edKWZPXs2Dh06hI0bN+LDDz/Eli1bcP/99+t1C6YmWPkCwOTJkz3q9OrVqz3OU/n6p7i4GAsXLsTOnTuxceNGOJ1OTJw4EefPn3elCdYvtLW14YYbbkBLSwu2b9+Ot956C2+++SaefPJJI27JVMgpXwC47777POrws88+6zpnuvJlhGKuuuoqtnDhQtf7trY2lpOTw5YvX25grqzJU089xYYNGyZ5rra2lsXHx7N3333Xdeybb75hANiOHTt0yqG1AcDef/991/v29naWlZXFfv/737uO1dbWssTERLZ69WrGGGOHDx9mANiXX37pSvPRRx8xm83GTp06pVverYB3+TLG2Ny5c9m0adP8fobKVxk1NTUMACsuLmaMyesXioqKWExMDKuurnalefnll1laWhprbm7W9wZMjnf5MsbY2LFj2aOPPur3M2YrX7LIKKSlpQV79uzBhAkTXMdiYmIwYcIE7Nixw8CcWZejR48iJycHBQUFmD17NsrLywEAe/bsgdPp9CjrgQMHwm63U1mHSGlpKaqrqz3KND09HSNHjnSV6Y4dO5CRkYERI0a40kyYMAExMTHYtWuX7nm2Ips3b0bPnj0xYMAAPPDAA/j2229d56h8lVFXVwcA6Nq1KwB5/cKOHTswZMgQZGZmutJMmjQJ9fX1OHTokI65Nz/e5Svwz3/+E927d8dll12GJUuWoKmpyXXObOUb8ZtGqs3Zs2fR1tbm8QMCQGZmJo4cOWJQrqzLyJEj8eabb2LAgAGoqqrCsmXL8P3vfx8HDx5EdXU1EhISkJGR4fGZzMxMVFdXG5NhiyOUm1T9Fc5VV1ejZ8+eHufj4uLQtWtXKncZTJ48GdOnT0efPn1w/Phx/PznP8eUKVOwY8cOxMbGUvkqoL29HY899hhGjRqFyy67DABk9QvV1dWSdVw4R3CkyhcAZs2ahd69eyMnJwf79+/HT3/6U5SUlGDNmjUAzFe+JGQIQ5kyZYrr/6FDh2LkyJHo3bs3/vWvfyE5OdnAnBFEaNxxxx2u/4cMGYKhQ4fikksuwebNm3HdddcZmDPrsXDhQhw8eNDDb45QD3/lK/bXGjJkCLKzs3Hdddfh+PHjuOSSS/TOZlBoakkh3bt3R2xsrI+H/OnTp5GVlWVQriKHjIwM9O/fH8eOHUNWVhZaWlpQW1vrkYbKOnSEcgtUf7Oysnwc11tbW3Hu3Dkq9xAoKChA9+7dcezYMQBUvnJ56KGH8OGHH+Kzzz5Dbm6u67icfiErK0uyjgvnCP/lK8XIkSMBwKMOm6l8ScgoJCEhAcOHD8emTZtcx9rb27Fp0yYUFhYamLPIoLGxEcePH0d2djaGDx+O+Ph4j7IuKSlBeXk5lXWI9OnTB1lZWR5lWl9fj127drnKtLCwELW1tdizZ48rzaeffor29nZXh0bIp6KiAt9++y2ys7MBUPkGgzGGhx56CO+//z4+/fRT9OnTx+O8nH6hsLAQBw4c8BCMGzduRFpaGi699FJ9bsSkBCtfKfbu3QsAHnXYVOWru3txBPD222+zxMRE9uabb7LDhw+z+++/n2VkZHh4cBPy+NGPfsQ2b97MSktL2bZt29iECRNY9+7dWU1NDWOMsQULFjC73c4+/fRTtnv3blZYWMgKCwsNzrW5aWhoYF9//TX7+uuvGQD23HPPsa+//pqdPHmSMcbYM888wzIyMtgHH3zA9u/fz6ZNm8b69OnDLly44LrG5MmT2RVXXMF27drFtm7dyvr168dmzpxp1C2ZikDl29DQwB5//HG2Y8cOVlpayj755BP2ve99j/Xr149dvHjRdQ0qX/888MADLD09nW3evJlVVVW5/pqamlxpgvULra2t7LLLLmMTJ05ke/fuZevXr2c9evRgS5YsMeKWTEWw8j127Bh7+umn2e7du1lpaSn74IMPWEFBARszZozrGmYrXxIyIfLnP/+Z2e12lpCQwK666iq2c+dOo7NkSW6//XaWnZ3NEhISWK9evdjtt9/Ojh075jp/4cIF9uCDD7IuXbqwlJQUduutt7KqqioDc2x+PvvsMwbA52/u3LmMMb4E+5e//CXLzMxkiYmJ7LrrrmMlJSUe1/j222/ZzJkzWWpqKktLS2N33303a2hoMOBuzEeg8m1qamITJ05kPXr0YPHx8ax3797svvvu8xnkUPn6R6psAbA33njDlUZOv1BWVsamTJnCkpOTWffu3dmPfvQj5nQ6db4b8xGsfMvLy9mYMWNY165dWWJiIuvbty/78Y9/zOrq6jyuY6bytTHGmH72H4IgCIIgCPUgHxmCIAiCICwLCRmCIAiCICwLCRmCIAiCICwLCRmCIAiCICwLCRmCIAiCICwLCRmCIAiCICwLCRmCIAiCICwLCRmCIAiCICwLCRmCIEzPvHnzcMsttxidDYIgTAgJGYIgCIIgLAsJGYIgTMN7772HIUOGIDk5Gd26dcOECRPw4x//GG+99RY++OAD2Gw22Gw2bN68GQDgcDhw2223ISMjA127dsW0adNQVlbmup5gyVm2bBl69OiBtLQ0LFiwAC0tLcbcIEEQqhNndAYIgiAAoKqqCjNnzsSzzz6LW2+9FQ0NDfj8889x1113oby8HPX19XjjjTcAAF27doXT6cSkSZNQWFiIzz//HHFxcfj1r3+NyZMnY//+/UhISAAAbNq0CUlJSdi8eTPKyspw9913o1u3bvjNb35j5O0SBKESJGQIgjAFVVVVaG1txfTp09G7d28AwJAhQwAAycnJaG5uRlZWliv9P/7xD7S3t+O1116DzWYDALzxxhvIyMjA5s2bMXHiRABAQkIC/va3vyElJQWDBw/G008/jR//+Mf41a9+hZgYMkoThNWhVkwQhCkYNmwYrrvuOgwZMgQ/+MEP8Oqrr+K7777zm37fvn04duwYOnfujNTUVKSmpqJr1664ePEijh8/7nHdlJQU1/vCwkI0NjbC4XBoej8EQegDWWQIgjAFsbGx2LhxI7Zv346PP/4Yf/7zn/GLX/wCu3btkkzf2NiI4cOH45///KfPuR49emidXYIgTAIJGYIgTIPNZsOoUaMwatQoPPnkk+jduzfef/99JCQkoK2tzSPt9773Pbzzzjvo2bMn0tLS/F5z3759uHDhApKTkwEAO3fuRGpqKvLy8jS9F4Ig9IGmlgiCMAW7du3Cb3/7W+zevRvl5eVYs2YNzpw5g0GDBiE/Px/79+9HSUkJzp49C6fTidmzZ6N79+6YNm0aPv/8c5SWlmLz5s145JFHUFFR4bpuS0sL7r33Xhw+fBhFRUV46qmn8NBDD5F/DEFECGSRIQjCFKSlpWHLli144YUXUF9fj969e+OPf/wjpkyZghEjRmDz5s0YMWIEGhsb8dlnn2HcuHHYsmULfvrTn2L69OloaGhAr169cN1113lYaK677jr069cPY8aMQXNzM2bOnImlS5cad6MEQaiKjTHGjM4EQRCEFsybNw+1tbVYu3at0VkhCEIjyLZKEARBEIRlISFDEARBEIRloaklgiAIgiAsC1lkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLCRkCIIgCIKwLP8fRAwazgqZP0EAAAAASUVORK5CYII=\n"
          },
          "metadata": {}
        }
      ],
      "source": [
        "file_names = []\n",
        "for i in tqdm(range(T_all)):\n",
        "  file_name = f\"fig_data_random/progress_{i:04d}.png\"\n",
        "  file_names.append(file_name)\n",
        "  plot_progress(y_best_data[:i+1], y_data[init_data_num-1:init_data_num+i], optim, file_name)"
      ]
    },
    {
      "cell_type": "code",
      "execution_count": null,
      "metadata": {
        "colab": {
          "base_uri": "https://localhost:8080/",
          "height": 488
        },
        "id": "ld6EjdnWLwWg",
        "outputId": "b7541c9e-9269-4cb7-e8eb-39e98ee19507"
      },
      "outputs": [],
      "source": [
        "APNG.from_files(file_names, delay=100).save(\"animation_random.png\")\n",
        "Image(\"animation_random.png\")"
      ]
    },
    {
      "cell_type": "markdown",
      "source": [
        "先ほどと違って一定のところで収束することがなくなったのでしょうか．実はこれは[大関研究室の森田くんが2023年に出版した論文](https://journals.jps.jp/doi/abs/10.7566/JPSJ.92.123801)の主張です．（[arxiv版](https://arxiv.org/abs/2309.02842)であればどなたでもご覧いただけます．）  \n"
      ],
      "metadata": {
        "id": "0Gwn3WLscAm0"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "本資料では$\\vec{a}$の推定に線形回帰を用いましたが，原論文ではもう少し複雑な推定方法を用いています．ここでは説明を省略しますが，具体的な手法を知りたい場合は原論文をご覧ください．日本語の解説としては[この記事](https://qiita.com/meltyyyyy/items/f92d911f551ceb32042a)や[森田くんが作成してくれた資料](https://github.com/mory22k/Note-2024/blob/main/202402/%5Bin%20Japanese%5D%20BOCS.ipynb)などがあります．"
      ],
      "metadata": {
        "id": "4CWY7AnPdREM"
      }
    },
    {
      "cell_type": "markdown",
      "source": [
        "また，線形回帰の部分をもう少し詳しく知りたい方は大関先生が過去に開催したQA4U第5回の[資料](https://gist.github.com/mohzeki222/23ac5ec8778edc082e58d1d2a060bc3d#file-bocs-ipynb)やT-QARDメンバーの鹿内くんが書いた[記事](https://qard.is.tohoku.ac.jp/T-Wave/?p=6252)をご覧ください．"
      ],
      "metadata": {
        "id": "J6Hl6bqwe8PV"
      }
    }
  ],
  "metadata": {
    "colab": {
      "provenance": [],
      "include_colab_link": true
    },
    "kernelspec": {
      "display_name": "Python 3",
      "name": "python3"
    },
    "language_info": {
      "name": "python"
    }
  },
  "nbformat": 4,
  "nbformat_minor": 0
}