{ "cells": [ { "cell_type": "markdown", "id": "3dae4474-c19d-4003-add1-821db4f66899", "metadata": {}, "source": [ "# 09-Finance Tutorial- Credit Risk Analysis\n", "### Modified for Quantum Rings toolkit for Qiskit 2.x" ] }, { "cell_type": "code", "execution_count": 1, "id": "7ab5133f-3fce-45d4-b4ba-7a71d373b53a", "metadata": {}, "outputs": [], "source": [ "# This example is from:\n", "# https://qiskit-community.github.io/qiskit-finance/tutorials/09_credit_risk_analysis.html\n", "\n", "# Modified for use with Quantum Rings toolkit" ] }, { "cell_type": "code", "execution_count": 2, "id": "2b7673a6-8f5a-45ac-a4b5-dce2465a68d0", "metadata": {}, "outputs": [], "source": [ "#\n", "# Setup your account\n", "# You can also save your account locally using the class method QrRuntimeService.save_account(...) and\n", "# invoke the QrRuntimeService class constructor without any arguments.\n", "#\n", "\n", "import os\n", "my_token = os.environ[\"QR_TOKEN\"]\n", "my_name = os.environ[\"QR_ACCOUNT\"]\n", "\n", "#\n", "# Set the backend of your choice, depending upon the task and your hardware configuration.\n", "# See SDK documentation for additional help.\n", "#\n", "\n", "my_backend = \"scarlet_quantum_rings\"" ] }, { "cell_type": "code", "execution_count": 3, "id": "0b6015ed-998d-448b-87c9-54723cc9f5bd", "metadata": {}, "outputs": [], "source": [ "import numpy as np\n", "import matplotlib.pyplot as plt\n", "\n", "from qiskit import QuantumRegister, QuantumCircuit\n", "from qiskit.circuit.library import IntegerComparator\n", "from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem\n", "#from qiskit_aer.primitives import Sampler\n", "\n", "from quantumrings.toolkit.qiskit import QrRuntimeService\n", "from quantumrings.toolkit.qiskit import QrSamplerV2 as Sampler\n", "\n", "#\n", "# Acquire Quantum Rings Backend\n", "#\n", "\n", "qr_services = QrRuntimeService(name = my_name, token = my_token)\n", "qr_backend = qr_services.backend(name = my_backend, precision = \"single\")" ] }, { "cell_type": "code", "execution_count": 4, "id": "a3e2b240-8f4b-49d9-8e3a-c369bfc0c367", "metadata": {}, "outputs": [], "source": [ "# set problem parameters\n", "n_z = 2\n", "z_max = 2\n", "z_values = np.linspace(-z_max, z_max, 2**n_z)\n", "p_zeros = [0.15, 0.25]\n", "rhos = [0.1, 0.05]\n", "lgd = [1, 2]\n", "K = len(p_zeros)\n", "alpha = 0.05" ] }, { "cell_type": "code", "execution_count": 5, "id": "0ec7de1b-3dc1-462e-8960-1d66cfb8ca03", "metadata": {}, "outputs": [], "source": [ "from qiskit_finance.circuit.library import GaussianConditionalIndependenceModel as GCI\n", "\n", "u = GCI(n_z, z_max, p_zeros, rhos)" ] }, { "cell_type": "code", "execution_count": 6, "id": "14185c8d-1ee2-4ca4-8c9f-d19fb1c7f742", "metadata": {}, "outputs": [ { "data": { "text/html": [ "
┌───────┐\n",
"q_0: ┤0 ├\n",
" │ │\n",
"q_1: ┤1 ├\n",
" │ P(X) │\n",
"q_2: ┤2 ├\n",
" │ │\n",
"q_3: ┤3 ├\n",
" └───────┘"
],
"text/plain": [
" ┌───────┐\n",
"q_0: ┤0 ├\n",
" │ │\n",
"q_1: ┤1 ├\n",
" │ P(X) │\n",
"q_2: ┤2 ├\n",
" │ │\n",
"q_3: ┤3 ├\n",
" └───────┘"
]
},
"execution_count": 6,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"u.draw()"
]
},
{
"cell_type": "code",
"execution_count": 7,
"id": "881a4f4f-b32f-4870-9a8a-a227eb164a5f",
"metadata": {},
"outputs": [],
"source": [
"u_measure = u.measure_all(inplace=False)\n",
"sampler = Sampler(backend = qr_backend)\n",
"job = sampler.run([u_measure])\n",
"\n",
"from qiskit.result import QuasiDistribution\n",
"result = job.result()\n",
"pub_result = result[0]\n",
"counts = pub_result.data.meas.get_counts()\n",
"shots = 0\n",
"for key, val in counts.items():\n",
" shots += val\n",
"quasi_dist = QuasiDistribution({outcome: freq / shots for outcome, freq in counts.items()})\n",
"\n",
"binary_probabilities = quasi_dist.binary_probabilities()"
]
},
{
"cell_type": "code",
"execution_count": 8,
"id": "1280c3d7-35db-44f9-940f-85c99f2e98d6",
"metadata": {},
"outputs": [],
"source": [
"# analyze uncertainty circuit and determine exact solutions\n",
"p_z = np.zeros(2**n_z)\n",
"p_default = np.zeros(K)\n",
"values = []\n",
"probabilities = []\n",
"num_qubits = u.num_qubits\n",
"\n",
"for i, prob in binary_probabilities.items():\n",
" # extract value of Z and corresponding probability\n",
" i_normal = int(i[-n_z:], 2)\n",
" p_z[i_normal] += prob\n",
"\n",
" # determine overall default probability for k\n",
" loss = 0\n",
" for k in range(K):\n",
" if i[K - k - 1] == \"1\":\n",
" p_default[k] += prob\n",
" loss += lgd[k]\n",
"\n",
" values += [loss]\n",
" probabilities += [prob]\n",
"\n",
"\n",
"values = np.array(values)\n",
"probabilities = np.array(probabilities)\n",
"\n",
"expected_loss = np.dot(values, probabilities)\n",
"losses = np.sort(np.unique(values))\n",
"pdf = np.zeros(len(losses))\n",
"for i, v in enumerate(losses):\n",
" pdf[i] += sum(probabilities[values == v])\n",
"cdf = np.cumsum(pdf)\n",
"\n",
"i_var = np.argmax(cdf >= 1 - alpha)\n",
"exact_var = losses[i_var]\n",
"exact_cvar = np.dot(pdf[(i_var + 1) :], losses[(i_var + 1) :]) / sum(pdf[(i_var + 1) :])"
]
},
{
"cell_type": "code",
"execution_count": 9,
"id": "08669ae5-c7d8-4684-b74a-3bb1e912af9b",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Expected Loss E[L]: 0.6592\n",
"Value at Risk VaR[L]: 2.0000\n",
"P[L <= VaR[L]]: 0.9512\n",
"Conditional Value at Risk CVaR[L]: 3.0000\n"
]
}
],
"source": [
"print(\"Expected Loss E[L]: %.4f\" % expected_loss)\n",
"print(\"Value at Risk VaR[L]: %.4f\" % exact_var)\n",
"print(\"P[L <= VaR[L]]: %.4f\" % cdf[exact_var])\n",
"print(\"Conditional Value at Risk CVaR[L]: %.4f\" % exact_cvar)"
]
},
{
"cell_type": "code",
"execution_count": 10,
"id": "51733958-b349-4600-b512-5ac8331e9e22",
"metadata": {},
"outputs": [
{
"data": {
"image/png": "iVBORw0KGgoAAAANSUhEUgAAAk4AAAHbCAYAAAA9GhWYAAAAOnRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjEwLjMsIGh0dHBzOi8vbWF0cGxvdGxpYi5vcmcvZiW1igAAAAlwSFlzAAAPYQAAD2EBqD+naQAAauFJREFUeJzt3Qd4E/UbB/C3g7aUvWVX9vqDbBCQoQwZMgRlqFQUVEBREASVoaCgIAoCDhRwgKIgyFQBAWXKkL0EKRvZoy20pc3/+f7SC2mbtEl6Gb18P89z5kgul7vc2b79jfcNMJlMJiEiIiKiDAVmvAkRERERAQMnIiIiIgcxcCIiIiJyEAMnIiIiIgcxcCIiIiJyEAMnIiIiIgcxcCIiIiJyEAMnIiIiIgcxcCIiIiJyEAMnIsqy1q1bJwEBAWrBelYQGRmpjjciIsLm69r5jBkzRrKSZs2aqePGI5GRMXAicuMv86z2y8+d8F1o34v1EhoaKoULF5by5ctL27ZtZdSoUfLHH394+3CJiGxi4EREXhUfHy8XL16Uo0ePysqVK2Xs2LHStGlTqVKlivz4449eOy4Gvxm3jhH5o2BvHwAR+Z9Zs2ZJ3bp11TrqjF+/fl0FT9u3b5dly5bJnj175ODBg/LYY4/JM888I59//rkEBqb9Ow/dQlmtTvmcOXPUYjRZpauUKLMYOBGRx917771SrVq1NM937txZ3nnnHVm6dKn06dNHLl26JF9++aXkz59f3n//fa8cKxGRNXbVEZHP6dChg2zatEly5cql/j1x4kTZuXOntw+LiIiBE5GvjvuZMWOGNG/eXAoVKiQhISFyzz33qMHT3377rSQlJaX7/iNHjsiLL76oWnUQfOD9xYoVk/vuu0+15MyfP1/i4uLSvC8xMVF1I7Vu3Vp9Ht6XJ08eNXD7wQcflHfffVcOHDggnoDPnDBhguXf1uvOzKpz9rvAeB7sT/PWW2+lGdCOsT/2jgHXBl2RuHZFihRRXYzW2zs7bmj16tXyyCOPSNGiRSUsLEzKlCkjAwcOlDNnzjg0ED899r4/7f1fffWV+veJEydsDux3ZVbdhg0b5Mknn1Tnj/PJmzev1KxZU958803VXevMsf7www/qvsT/I9mzZ5eKFSvKsGHD5MqVK+keA1GmmIhIV2vXrsWgG7WMHj3a6fcfP37cVKlSJcs+bC2NGzc2Xb582eb7f/jhB1NISEi678eyd+/eFO+7efOmqUmTJhm+79FHH3Xpe8F3oe0D35EjYmJiTHnz5lXvCQ8PN8XHx9v9rm3t05XvonTp0hlu37t3b5vHsHLlStNDDz2U7vZYx3P4HFus750xY8bYPYY8efKY/vjjjwy/6/TY+/6s35/eYq1p06bqOTzakpiYaBowYEC6+8M5/fbbbxke65o1a0xPPPGE3f2UK1fOdO7cuXTPnchVHONE5EOio6PVX9D//vuv+nenTp1UqwhaSI4fPy7Tpk2T9evXq7/a0Z2FaftBQUGW9//333/y9NNPqxYrTPFHy0SDBg2kYMGCcuvWLTVzDe9fvHhxms9GK8Off/6p1tu3by+9evWSUqVKqVaBCxcuyN9//60GbmfUiqGn8PBwuf/++2XFihUSGxuruuvq16/v0Htd/S5+++039Z7//e9/6t8vvPCC9O/fP8U2+fLls/mZr732mhrYjhYitCyVLl1aHceNGzecPvfly5erwfJaK0r16tXVIHrMNJw5c6Zax3Xat2+flCxZUvSE8+3atatqBfr555/V/ffrr79map/Dhw+X6dOnW8a44buqVauWxMTEyJIlS9S9rZ3TX3/9JTVq1LC7r5EjR6quXPz/8dRTT1m+Z+wf3xuu7SuvvCLfffddpo6ZyCaXQy4i0r3F6dVXX7W8980330zzelJSkqlXr16WbWbMmJHi9S+//NJui5K12NhYtVgrWbKkel/Xrl3TPUZ7LV3uaHECfA/a+77++muHW5wy812Ao9fQ+hjsXTdrjrY4YalVq5ZqCUwN34O2Tbdu3XRvcXL0WB1tcdqzZ48pMDBQvV6tWjXT1atX02yD1jptm3r16qV7rFjGjRtn8/+PVq1aqdeDg4NNFy5cyPC4iZzFMU5EPgLjbL744gu1XrVqVZv5g9Dag7FPBQoUUP/GX+nWzp8/b2kRsTVrTYPxIFhsvbdJkybpHidmuHmSdq5w9epVh9+Xme/CVRUqVNA17xPSMOTMmTPN8xgj9PDDD6v1RYsWWc7VV33yySeWcXm4xzGuKbU2bdqo1lVAi9O2bdvs7q927dry+uuv2/z/Y/DgwWr9zp07snnzZh3PgsiMgRORj9ixY4dcu3ZNraObx7oLzlru3LlVfiPAQO1z585ZXsMAYi3AQBeLM7T3YrA0usV8hXXgcPPmTYffl5nvwlWPP/643evmLHQVIkCwRwsyECD4eg4lDHDX/iBIr6u1b9++ad5jS8+ePe12GVt/Z1qXN5GeGDgR+QiMVdFkNI7H+nXr92FsjfbXPHIitWjRQj788EMVlGHGXHp69+6tHjF2BGNQMCYIrRnpzXTyBOtgCUGjozLzXbgK45D0oiUItadevXqW9b1794ovt6T+888/Dt3XmF2XLVu2NPd1apUqVXKoRdSZQJvIUQyciHyE9RRqDGZOD1IF2HofurUw0LZ48eIqo/batWtV10WdOnXUL5QuXbqoAd72BtyiFQN/yWMwOAbaYnscC7q6Ro8erQbgehqSYLrSTZiZ78JV9gaNuyKjewCpDjS+PP3euns1o3NC0KR1zaZ3Tpg0YI91hnl3Bcjk3xg4EfmgzMxcwxglzCpCvid0aZQoUUI9j5ldaEHCbDyMJ0ndHYdfWsjSjb/0MZsKs9mQ8wj2798vb7/9tpQrV85j3V4azObTYIaZJ74LV+nVTQeenL3oKUY8J/I/DJyIfIR1a0pGLTvWg4FttcIghQDSCcydO1dOnTqlxnp8/PHHavAyYGr5G2+8YXPfKK6LQrsbN25U08NXrVqlpvUjKEC6hB49eqQYV+VOCGjQdaiNdULSSmdl5rvwpozuAevXU98D1q0u6SVLRSoAd7NuhcvonDBe6/Lly16ZhEDkKAZORD7CeubX1q1b090Ws45svc8ebcwSZipprS7IuuxI0PHQQw+pTNgoewLIgaR3F5c9s2fPVsEbIL9PcHDmU8+5+l14WnqzylK/nvoe0ErVZDQTEVnV3d1CFBoaqrLAO3Jfo3UxISHB4fuayBsYOBH5CMwG0gYzo9SFvZYCDHjVftGjdUibPeYIDK7WBh1bjx1yBBJzapx9ryswoHjEiBEpEijqKaPvAkEj2CpN4wkY8G3dTZkagllAS2DqMicIDjVIomnP999/n+4x6PUdIPjWunytg/7UtHQc1u8h8jUMnIh8BP4yf/bZZ9U6xhmhuyw1DHJGa4n2ix7r1tDtlF43GlpvtF9c1r9cMRB36dKlav/2IKO2xvq97oAWLYyx0mZFIYBKL5O0La5+FxotID127Jh4S79+/Wx2p82bN09lUwdkz04dPOO701rnMJPQ1nVFC2J6QQxo+8VkgczMUEP2da37EOdkK5M67i+MsdNmDGY0q5DIW1hyhciNdu3apYrmZgRT5VHeZNSoUfLTTz+pcThIpIhWB4wvwi8wreSKlrOnYcOG6peQNZSYwIDnli1bSqtWrVR3B8aK4JcegjG8XysO+/zzz1veh19kmL6PwquYbYZp4yhjgV++CD4QVGmtAZilhm6zzMC5oPQJ4Jc6Ph9pD9A6gs9C2RINzvGdd95x+jNc/S6sgw8cJ2bmffbZZ9KoUSNLCwxaqzKaIZZZmP2H7wOPKE+CvE4I9hYsWKCOR+uSmzRpUpr34ti6deumvgMEkLi2AwYMUDPxTp48Kd98840sXLhQnaM2hswWvA5o/cR3hGLJ2nUDTBZwBI59yJAhKljbvXu3KrWCc0L6AQSGuOZTp05Vs+AwIUE7PyKf5HSucSJKV+rSEI4sixYtcqrIb6NGjWyWPtFKZGS0PP/886roqvVnOvK+okWLmrZv3+7S9+Jo4VhtqVKlimnhwoUulwxx9bvQ/P3336bQ0FCnivw6UkrGmSK/6X1nuXPnNq1bt87u55w/f95Uvnx5u+/v3r27afXq1ekeO76XBg0a2N2Hs0V++/fvn2GR319//dXm+535nh0tl0PkCrY4EfkYtPrgr3IUckVBV7SOoEUGrSX4Cx0zxDC13nrmlAbdMmhh+f3331WrDVqL0JKDcTAoBItWKnQHNm7cOMX70LqEbht0/6AF4sSJE2oGFGbRYdwVxlKh9QatP84koXQE0iBgn3ny5FHpBtDCghai1MfoLFe/Cw1m8KFkB1pJMMMQ34enxzuh1RHHiVmAaH3CQG8U3G3btq3qvtQGt9uC1iUMxn7vvfdUKyZamnLkyKFa3nAdcR9llHEc9xi60N5//33VKoRuS7QQpdelm96+kBuse/fuqkUJBaXxnaKLukyZMuqcXn75ZSlUqJDT+ybypABETx79RCIiIqIsioPDiYiIiBzEwImIiIjIQQyciIiIiBzEwImIiIjIQQyciIiIiBzEwImIiIjIQczj5CJk0j179qzK3KtHIUwiIiJyP2RhQgUB5ESzlQ8vIwycXISgCUn0iIiIKOs5depUuklk7WHg5CK0NGlfvN6ZlMm/xMTHSLEPiqn1s0POSo6QHN4+JPIXd2JEfjLfe9LlrEgw7z3DQHHoYsnX9uxZkRy8thpUYkDDh/Z73FkMnFykdc8haGLgRJkRFB8kYq4dq+4lBk7kMXeCRMKT1/FzjIGTcQQF3V3HtWXglIarw2w4OJyIiIjIQQyciIiIiBzErjoiLwsODJbeNXpb1ok8JiBY5N7ed9fJOIKDRXr3vrtOugkwYV4euTS4LE+ePHL9+nWOcSIiIvKT39/sqiMiIiJyENvviLwMjb6xCbFqPTxbOBOqkuegwyHRfO9JUDimGXn7iEjPaxubfG3DeW31xBYnIi9D0JRzfE61aAEUkUcgaPohp3nRAigyBgRNOXOaFy2AIl2wxYmIiNwqISFBEhMTvX0Y/iUuTqR06bvr1nmdDCQoKEiyZcvm0c9k4ERERG4bhHvp0iWJwy9u8qykJJFPPzWvnzsn4kJNtqwiNDRUChYs6LGJWgyciIjILUHTmTNnJGfOnOqXGloFOH7Pg9DCd+uWeT0iwpAtTiaTSbVmYnYc7jXwRPDEwImIiHSHliYETSiiyoDJC6y7RsPCDBk4Qfbs2VXNudOnT6t7zhOBk3Hb7oiIyCvQCoDuOeTKYdBE7oZ7DPca7jnce+7GwImIiHSlDQT39KBd8l/Zku81T0xCYFcdkZcFBQZJ1ypdLetEHhMQJFKy6911vXfP1ibvwXefL9/ddYML8OA5MnAi8rKw4DD5sduP3j4M8kdBYSJNeO8ZEmbRlS3r7aMwJAZOPihi+HJvH4LfiprQztuHQEREPoxjnIiIiDzUnZTREhkZadk+Kioqzevbt29Psc9mzZqp59etW5fh51+7di3N/hx5H6XEFiciL4uJj1HlViB6RLTkCMnh7UMif3EnxlxuBR6LFgnmvecJvXv3tvta48aN0zxXpEgRadOmjVpHTiyHYJD033+b12vWVOkIQkJCLJ+9YcMGOXbsmEvH7+8YOBEREXnQnDlznNq+UqVKTr/HlvDwcMt+0LLFwMk17KojIiIichADJyIiIiIHMXAiIiIichADJyIiIiIHcXA4ERF5ZTapPcigj8SwjmwbGBAo2bNld2nb2IRYMZlMNrfFVP3wbOHi6SzXixYtkk6dOrnlc0kfDJyIvAy/JNqWb2tZJ/IYlFkp1tZtJVfSo6XgsAX/PyzveTcRcOFJhVWQY0vT0k1lXeTdXEQRUyLkUuwlm9vWKVZHtvXdZvl3lelV5MT1Eza3rVKoiuzvv188nY6gVKlS+nwIgrM8ee6uk24YOBF5Gf6ytv4lQeTRkivNeO95mh6pBRwquVK+vPs/xw8xcCIiIo9Dsld7Ure8Xnj1Qrrdb9aiBkU5vO2BAQfS7aojsoWBExEReZwzGfLdta27xjCRsXFWHZGXYTBrjndzqCW9ga1Ebim5Mj+HecE6GQdKruzcaV6wTv4RON26dUtGjRolFSpUkLCwMClWrJj06dNHzpw549L+UDDx+eefl3vvvVdCQ0NVzZ+GDRvKxIkTdT92Imdg4Ku9wa9EbpUYa17IeJKSzAv5R1fd7du3pUWLFrJlyxYpWrSodOzYUQU+s2fPlmXLlqnny5Qp4/D+Vq5cKV27dlXBWK1ataRBgwZy+fJl2bt3r3z22WcydOhQt54PERGRVicuvVl1b7/9ttP77N+/v+TOnfvuExi7FWsOiGs1biwzPvnEtYOlrBM4jRs3TgVHaBH67bffJGdO89TVyZMny5AhQ1TL07p1d6egpufQoUPSpUsXyZUrl6xatUruv/9+y2tJSUmyE02ZREREHvDVV1/Zfa1GjRouBU4HDx60+1pY/vxO74+yWOAUHx8v06ZNU+vTp0+3BE0wePBgddOtX79eduzYIbVr185wf3gPWrAWLlyYImiCwMBAqVOnjhvOgoiI6C57M/gyw24DAsY1/f23eb1mTd0/15/55BinjRs3yvXr16Vs2bJS08YFR5cbLF26NMN9nTp1Sn799VfVrde2bXKiNyIioiwCvSbo3sOCISuuio2Ntexnw4YNuh6jP/HJFqfdu3erR4xFskV7fs+ePQ5F4+iOQ0vTnTt35KefflKBWWJiolSrVk0ef/xxyZcvn85nQEREpI///vvP0r03cOBAiYiIcLk3J71uQsrCgdPJkyfVY4kSJWy+rj1/4oTtVPnWDhw4oB7R3dekSRM1bsraG2+8IQsWLJDmzZvrcOREzkNSPpSN0NaJPCdQpHBTX+6A8GsIkFzu3kMCz1y57q6LSN68ed3SXehvfDJwio42Z5QND7ednCxHDnOCs5s3b2a4r6tXr6rHL774QgVP8+bNkzZt2sjFixdl7Nix8u2330rnzp1l//79Urx4cbv7iYuLU4vmxo0bTp8XkS0oOmpda4vIY4KzizzEe8+QUHKlYkVvH4UhGf5PDHTTAbrpkHagR48eqmsOuaG++eYbqVu3rhpPNWPGjHT3M378eMmTJ49lKVmypIfOgIiIiHyFTwZO2iw6DGSzJSbGnOEW6QUc3Rceu3Xrlub1p59+Wj1ill56RowYoQIsbcGgcyIiIvIvPtlVhwRgcPr0aZuva8+XLl06w31p22Cftoo2aoPsLlywX0QSkGkcC5HeUGYlYkqEpUCpM7W2iDIFZVZ+Th5o3DFKJJj3nmEgHcHeveb1//1PJChl4WQyWOCEBGBgLzGl9nz16tUz3JeWzkAb65TalStX1KN1rigiT7sUe8nbh0D+Ko73nmHduePtIzAkn+yqa9SokRpHdOzYMdm1a1ea1zELDjp06JDhvpCGoECBAnL+/Hk5fPhwmte1Ljpb+aKIiIiIfD5wCgkJUbkqYMCAAZYxTVrJFeRvatq0aYqs4cg0XqlSJTUWyVpwcLDKHI4pmNiX9Wy41atXy5w5c1QX3nPPPeeRcyMiIqKsyye76uDNN99Ugc2mTZukfPnyKgcT8jZt3bpVChUqJLNmzUqx/aVLl1SL0rlz59LsCwV8165dq/aH2XQo8IvtkdMJiTDfeecdqVevngfPjoiIiLIin2xxgrCwMBXsjBw5UuVzWrx4sQqckCoeY5xQQsVR2bJlkxUrVsh7770nBQsWVCVY9u7dq1qtULbl9ddfd+u5EBERkTEEmJhG1CXo8sM4LKQmyJ07t677jhi+XNf9keOiJrTzyqy6nOPNkxOiR0RzVh15dlbdD8kTYx6L1m1WHYqqHz9+XO699171RzB5QeoivwafVXfbiXsus7+/fbbFichfoMxKnWJ11MKSK+RZgSL565gX/jpwm549e6qxtKhWkZG//vpLbVukSBGVuNkZ6JHBe9USHCwBdetK9saNpXylSmocLwKLjLRo0UKVNbOulIHCwtinozXyXn75ZcmePbulfJrR8P8UIh8oubKt7za1YJ3IoyVX2mwzL1gnt3jyySfV49y5czPcFmXAAFUuMLnJ1ZnpvXv3VsuDDz0k165dk88//1zuu+8+NbnKnuXLl6shMhi+kpm8ha+99pplrLIRMXAiIiJyo1atWqkWJExg2rZtm93t0MI0f/78FMGWK5599lk1YxzLsmXL5OjRo2pML7qoMMvcHgRMmHyF92dG0aJFVdCGIPDAgQNiNAyciIiI3CgoKEi1IFm3KNny22+/qSoWlStXTpFuJ7MwngeTo7TchRgPlNrGjRtVa9Tjjz+uUgJl1hNPPKHSAH366adiNAyciLwsNiFWIj6KUAvWiTzmTqy55AoWrJPbIJAAtCghDY4tWlcetkX32scffyytW7dWpcPQdYZkzm3atJFVq1Zl/IH4DHTLYUlMlKpVq1patWxV0vjiiy/UoxbgZVajRo1UqTMEirYCtayMgRORl+GvshPXT6iFk1zJs0wiMSfMC9bJbdCChJak//77z2bgg0TPP//8sxqE3atXL5Vn8KWXXpIjR45IxYoVpXPnzuoRrVIIplLnMrQpPt68iMjNmzfVY2BgoArAbI1vwoBuvXIaBgQEqO5BBGnIx2gkDJyIiIg8QBu3ZKu77qefflLB0wMPPKBamBAkbd68Wc2EQ7D0/fffqwBkx44dquvtlVdekejoaIc/+5dffrHMmkvdFXfo0CG5ePGiKj3m6oB0W7QgTCttZhQMnIiIyDs5pOwtibcd3/bOrUxsG5vOtvp3XaIlCS0xSOhsXUrMOpjSuvSQjwhVLlJDcKOVD8MMuIxcunZN5s6bJ6+++qoa+D1lypQ022gz7RCs6alSpUrq0VbN2azMZ0uuEBGRgWmJN20p1lakmVUi4IWFRRLtBDKFm4o8tO7uvzFeK+6S7W2RrwqpFzTLqyR3U9qQp4pIu/2iJ4z5QYsSWmAQPCGQAnTfrVmzRiVu7Natm2V7jIXC82hpQjkxLbfSP//8k+Ixtaefflot1tCKhQHgxYsXT7M9BqRDvnz5dDxbkfz586tHtGYZCQMnIiIiD3bXIXBCC5MWOH333XcqSOrSpYvqhoPTp09L+/btZffu3Xb3pY1bsjUwu1zZspJ06ZKcuXBB1v/9typZhhQBKDmGWX7WkEEbcuXKpeOZiiUrNwa6GwkDJyIi8jyUeLEnIFV5kEfNLSIOjTjpGOX4tu0OpDMoPkDcoWvXrjJw4EBVdB4tPYULF7Z001nnbkIuJQRNjz76qAwbNkx1oyGwweBuJLNEJnB7k0nw3kjsK7nkyoHQUGnaooVqvfrwww9Vt501LVizF4i56npyQJY3b14xEo5xIvIyjHmoUqiKWrBO5DkB5i4pLG4KFOxCXTx7S1CY49umznju1Lbh6Wwb7pbTRpDyyCOPqLQAaGnCwGwM+EYBeqQaAIx/wsw7JM1E+gIMssb7EDTBv//+69iHoWZbWJhUqVJFpk6dqp569913LQGNBsEbXLlyRddzvZqc9gBjq4yEgRORl4VnC5f9/ferBetEHoPgAON4sLgpUKC0tAHgyNuk5W5C4sls2bKpdQQ2SUlJKgN36m61hIQEWbRoUcYfgvdVq2ZegoKke/fuquQKgpnp06en2LRGjRrqEZnN9XTw4EH1iM81EgZOREREHoSWJbQwofyKllnbupsOLUBoYdq3b58a0K3BOCjUgUNuJ2ehNXvMmDFq/aOPPpLY2LuD7dENiM/E7DdnCwtnVLAYkM/JSBg4EREReRBaltACBJcuXZLy5ctL/fr1La8jlxLGNSGIQdCBWnfYvly5cirQQjoCV3Ts2FFq1aqlZrnNnDkzxWtt27aVW7duydatW+2+HzP7GjRoYHdBEk0Nxl9hEDzGN91///1iJAyciLwMZVaqzqiqFpZcIY9CrqLlVc0LS654lHULk9Z1l7rg7ldffSXVq1dXrU4YTI4uNWQUr1OnjmMlV/btMy9WJV60VqdJkyZJfHJWcejbt696nDdvnt1dYvutW7faXazTDmzYsEFOnTqlzhNpFowkwMQaDy5B8jE0paIvWptyqZeI4Vb5S8ijoia08/hnxsTHSM7x5pw20SOiJUdIDo8fA/kpJHrU8ilhlhsGResAtcmQ8RpJHI32SzPLQLCUPKtOatY0j3nKAJJrIg0CFtTGy4znnntOtWrt3bvXUifPnZy55zL7+5stTkRERCTvvPOO6jpM3Y3nrHPnzsnXX3+tWtI8ETR5GgMnIiIiUuOcmjdvLhMmTLBkKXfFe++9px7HjRsnRsTAiYiIiJTff/890111H330kRpojhIzRsTAiYiIiMhBDJyIiIiIHMRadURehsR0pfOUtqwTeU6ASA7zvefxkivkfiEh3j4CQ2LgRORlKLMS9XJ6hUmJ3ARlVtItiktZFtIPVK/u7aMwJHbVERERETmIgRMRERGRgxg4EXnZrYRbUndmXbVgnchj7twS+aWuecE6GUdSksiBA+YF66QbjnEi8rIkU5JsP7vdsk7kOUkiV7bfXSfjQDW12OT6g6yspiu2OBERERE5iIETERERkYMYOBERERE5iIETERGRB8XExMjkyZNVQd0iRYpISEiI5MuXTxo2bCijRo2SkydPypEjR1RC3Fy5ckmsNlYpgwK92P7jjz92+niioqLUe62XoKAgyZ8/vzRt2lTmzJkjpgzGSf3xxx/qfdOnT0/xfGRkpHoe+8jIuXPnJHv27NK/f3/xZRwcTkRE5CGbNm2SRx99VM6fPy/h4eHSoEEDFTxdv35dtm3bJlu2bJH3339fli1bJvXq1ZO//vpLfv75Z+nRo4fdfV64cEFWrVolwcHB0r17d5ePLUeOHNK1a1e1npCQIP/8848KiLCsW7fObvBjMpnk1VdflRIlSsizzz7r8ucXLVpU+vXrJzNmzJCXX35ZKlSoIL6ILU5EPqBgeEG1EHlcaEHzQm63a9cuefDBB1XQ9Nprr6mAZ82aNTJv3jxZvny5en7hwoUqADl9+rQ8+eST6n3ffvttuvv9/vvv5c6dO9KmTRspVKjQ3ReCg82LgwoWLKiCIyxz585VQduiRYvUa1999ZVs2LDB5vsWL16sgr7BgwdLaGioZMawYcMkKSlJRo4cKb6KgRORl+UIySEXh15UC9aJPCY4h8ijF80L1slt0CqDQOj27dsyZswYmTBhgmrhsRYYGChdunSRHTt2SJ06dVTrUbZs2eS3336Tixcv2t23Flg98cQTKUuu3HefecG6izp16qQCMvj1119tboMWInTt9ezZUzKrePHiqgsTAdt///0nvoiBExERkZv98ssvsm/fPtWa9MYbb6S7bZ48eaRatWqqBah169aqNWn+/Pk2t0V3Glp7cufOLY888oh67s8//5SBAwdK9erV1dgpjBuqVKmSDB8+XK5du+b0sVetWlU9ooUstePHj6tWsxYtWqguRz0gAENXoSPjoryBgRMREZGboSsOunXrpsYiOSqj7jrteYxNQoAEQ4cOlS+//FL9G12DWG7cuCHvvfeeNG7cWKKjo5069ps3b6rHwoULp3ltxYoVqjWtWbNmohdtX9p35ms4OJzIy1Bm5eG5D6v1lb1WSvZs5h9+RG6HMivrzPeeNFspEuzBey8mxv5r6FoKC3Ns28BAkeSAweltMVvN3myxgACR8HDRc3wT1KpVy6n3oRUJLVBbt26Vo0ePSrly5VK8jrFI1gEWjB49Wu5v0EDyaC1E5ctLXEKCvPTSS/L555+rGX2YvecItPygRQm0LjtraN2CunXril7KlCmjWtswxgpdm2HW94IP8OkWp1u3bqmLi5H1+OKKFSsmffr0kTNnzji1n4iIiDRTLa2XQ4cOue0ciDKCMivrT6xXC0uukGcliVxYb148XXIlZ077y6OPptwWLR32tn04OfDTRETY3/aBB1JuW6WK/W11DATg8uXL6jHF4G0H4HefNtMtdavT5s2b5dixY1KyZEmVNkDz8MMPS57cudFUZF5MJjVo+6OPPlKtXZil50jAdODAATXOCp8xYMAAadSoUZrt9uzZox4rVqwoesL+4uLi5ODBg+JrfLbFCVEm+kwxNRNTFDt27KhyTcyePVtN08TziEqd0bt3b5vPI5onIiLyRRj0ja43tC5hYHnq1qZevXqpRgBraGBYunChHIqKkhthYZKU3LKGnFEYF2XLiRMn0uwHxo0bZ3dc1oXkVi2MpdITckhBeoPivcVnAydcKARHSAiGGQU58ReAiGpiHDJkiGp5Ql4JZ/jqQDMiIr+T3jib1LPAbAxKTtH9Zi0qyvFtDxxIv6tORwUKFHA5EEBrUqlSpVRXHbrs6tevn2LAuHU3nfZ7EgPB0Wrkah4nJOnEoHMEU2+99ZbKKdWyZcs070H+KdB+R+sFg93BlcHsftlVFx8fL9OmTVPryEJqfUGQJwIzBdavX6+mbBIRURaEqfj2ltRjWtLb1nrMkrPbYgyTvW11HN8E9yEtgIjs3LnT6feiFQitStbddZild+nSJTVmqgq6HJOhwQGNC0iuOWf0aIlaskRux8SoAdxY0IPjSB6nH3/8UXXRvfjiiyoAe+qppyyDxG312Dg74DwjWkCWN29e8TU+GTht3LhRfWlly5aVmjVrpnldi4iXLl3qhaMjIiJyTrt27dQjAhK0FjlLa1VCKxPebzN3k4glYeU7Y8dK7/btpXTRopaklBg3jCSbjkJuJrReIR0B3vfhhx+m2UabaXflyhXR09WrV10aE+a3gdPu3bvTnX2gPa8NSnPUxIkT5fnnn5dBgwapmQW+2HdKRETGgxlpCECQEfydd95Jd1ukDti/f3+K5ypXrqx+9+H3FrKLL1myRAU2qUuxaAEH8kWlhqAto5pzqWEwOYbOwJQpU9K0LNWoUUM9Hj58WPSESVsI+HDevsYnAycUOLR34a2fR9+rs6ncP/vsM5k6dao899xzarbdrFmzdDhioswJzxauFiKPCwo3L+RW6G5DKxFmyWGA94gRI9Q4ImsIahAQIWs4xhfZa3VCcku0HmHM0T333JNiG62+25ezZklCYqJlXBdmyKHMiyswOQu9P2hV+uSTT1K81qRJE/Vo63hdhS5CzELEuCpfS0Xgs4GTFtGij9YWLU29rf5We3kwfvrpJxVooco0srdirBSmOqIgoSNTM7Et/gqwXoj0gDIrMa/HqIUlV8ijUGbl8RjzwpIrbodxTqtXr1YZtlFyBd1cDz30kBq/1L59e8sM8lOnTqkUA6mhdQmtTBjbZGtQODz99NMqmFq6bJlU7NFDHn/vPWnZpo36bAQ5pUuXdino02bzTZ48Wc16t059gNczmqw1duxYVdDY1tK5c+cU22r70ro3fY1PBk56QwsTLgxmJSCTKppLP/jgAxU5I8J3JAofP368GgSnLbZuaiIiovQgFxJmx02aNEkljcSQkx9++EGN7UUvCJJXIl0Asn2nhoCrVatWah2TplBHztbsPbT+oGwJJlqhBQupCRC4fPfddy4fNxogateurcY6WffU3HvvvSr4W7t2bbrjp/799181I9DW8vfff6fYFkWPUaMvMjJSfFGAydkOTw9AaxAGob3yyisqurU1BgrRM/p7MzOzDhWYEeEjDwXq7eCmTa/FCYsGLU4InjCIXZs2qZeI4b6ZZt4fRE3wzb9wiLIStEjgZyp+qfpiVwvp6+eff1ZBHIJBzOjLDIwBQ6sYJoHZq8+X2XsOv7/RAOLq72+fbHFCy5D2BdqiPe9Kk2PqStSYuQfnzp1Ld1sMUsMXbL0Q6eH2ndvSbl47tWCdyGMSb4usa2desE7GkZSECsDmBetu1LFjRzUeCQ0e1g0MrsAkLvxufvvtt8VX+WTgpI3St5fvQnse+ZwyS5uBoI2bIvK0xKREWfHPCrVgnchjTIkiZ1eYF6yTcaAzCbmQsHigY2nixImqS3DmzJku7wMNGJjx3rdvX91LuBg+czj6gNGMhpH1KIyoJQ7TLFiwQD126NAhU5+D6Z6YQolB6JUqVcrUvoiIiPzVAw884HSqg9QwdAazBX2dT7Y4oZYOplsCCgtaT9nEmCcMpkMKegxU0yDTOIIfTPG0tmLFCvn999/TfAb20a1bN3WhMbMOn0lERESU5Vqc4M0331TTNjdt2iTly5dX0yiRTgAj8JFJNHX+JUzPROtR6rFKf/31l6qzg/FQ6AJE6xJG96O7D9lXmzVrpqaFEhEREWXJFifAqHhMbxw5cqQKdhYvXqwCJ0xPRNBTpkwZh/bTunVrVRAYg7kx3RPdfJgK2rhxY9UXi+AMKQqIiIiIsmyLEyCgwch6R0bXIzmXlqDLWsOGDdVCREREZNgWJyIiytp8ME0gGZTJg/eaT7c4EfkDlFkxjeYvGPIClFnpqf+9h7IgkJCQwKEQ3oJrUKeO+IuEhIQU9547scWJiIh0hXIZSBqMzMxsdSJ3wz2Gew33HO49d2OLExER6a5gwYIqISIqPSAvH36hoRgskZ4BE1qaEDRFR0dL8eLFxRMYOBF5GcqsPLnIXOX8m87fSFgwa3uRh6DMyibzvSf3fyMSpN+9p5WlQqoYBFDkYWjpu3TJvF6woIiBg9bQ0FAVNHmqFBoDJyIvQ5mVBQfM2fDndJzj7cMhf4IyK6fM956Y9L/3tLqeaBVITGRJF4+KjRVp29a8jjJl4eFiREFBQR7pntMtcIqPj5eDBw/KxYsX5dq1a5I3b16VnLJy5crMxE1ERAp+sXn6l5vfQ6B64oR5PTQUyRG9fUSG4XTghCBpzpw5snz5cpWV21YlZDSboVJy+/btpXfv3iqYIiIiIvKbwAnZtpHFe9GiRaqlSRv8h3px+fPnV82xGKB19epVOXTokPzxxx9qQemULl26qCSW5cqVc+e5EBEREXk/cELBXZQnQR918+bNpWfPnqrG27333mv3PagHh5Ip8+bNkx9++EEWLlwo/fr1k48//ljP4yciIiLyGIfyOKGg7gsvvCAnT56UVatWydNPP51u0ASoJffMM8/ImjVrVI25559/Pk1hXiIiIiLDtTih9eiee+5x+UMwTXDKlCkyYsQIl/dBRERElCUCp8wETe7YD5GRhGcLl+gR0ZZ1Io8JChd5LPruOhkH0g9EJ19bg6Yi8BbmcSLyMmRTRr06Io9DUkTUqyNjXtscvLY+Xatuz5498tRTT0mdOnVUKoI+ffqoHE9ERERERqFL4PTjjz+qtAQ///yzyuIZGxsrX331ldSoUUN++eUXPT6CyLDi7sRJ5OJItWCdyGMS40Q2R5oXrJNxIMdiZKR5sZFvkVwXYNKhdDVm2FWtWlW+//57yZkzp3ru77//lhYtWkhERIRaN5obN26owpXIXaV3fZyI4ct13R85LmpCO49/Zkx8jOQcb/7/BmOd2G1HHnMnRuQH872nxjqx2844YmJEkn8fq7FO7LbT7fe3Qy1OyOFkz+3bty3pBrSgCWrWrKkCJ3bXERERkVE4FDghKKpfv75s27YtzWthYWEqclu3bl2K52NiYlRLE2fSERERkV8FThs2bFDVrRs2bCh9+/aVS5cupXi9f//+MnnyZHnooYdk+PDh8tJLL6muu6ioKPUaERERkd8ETgiYduzYocqloFZdxYoVZcaMGaINjxo3bpxMmjRJdcu9//77Mm3aNElKSlKPw4YNc/c5EBEREfnWrDrkmkHZlSNHjkjXrl1VqxJm0m3atEm9NnjwYDlz5owabIUF5VnY2kRERER+nY4gf/788tlnn8mWLVskJCREmjRpIpGRkXLx4kX1eq5cudRCREREZDQu53FCoksET5hxh1xN5cuXV/Xo0EVHRI5DmZULr15QC0uukEehzEqXC+aFJVeMBWVWLlwwLyy54r3A6b///pPff/9dFi5cqGbYxcfHqwzhhw8flieffFJeffVVue+++2T9+vX6HiWRgaGru1COQmrBOpHH4H4LK2ReeO8ZC65noULmhdfW84FTXFycGq9UqlQpadmypXTr1k0aNGgg5cqVkwULFqh0BBg4jgHkefPmVfmbevbsKWfPntX3aImIiIh8PXAaOnSofPrpp9K8eXOZO3eurFy5UqUfCAwMlO7du8v27dvVdtWrV5c//vhDvv76a9XqVKlSJZk4caK7z4EoS0OZlQHLB6iFJVfIo1BmZdsA88KSK8aCMisDBpgXllzxfMmVwoULq9YmLUDS7N27V9WjGzJkSJoAKTo6WsaMGaNSEiC7uNGw5IoxseQK+RWWXDEullzxbskVZAEvUqRImue1rOC3bt1K8xrKryC3065du5w+KCIiIiJf5FDghC66X3/9VQVCFy5cUFnE9+/frwaGYzBr06ZN7b4X3XVEREREfhM4TZ8+XSpUqKCygBctWlTVp8N4phUrVqgSLBgsTkRERGR0wY5sVLp0adm3b59KQ7B79265evWqGvP08MMPqwCKiIiIyB84FDgBZtChZYmtS0REROSvXM4cTkRERORvHAqckLdJDxgTRUQpZc+WXY4POq4WrBN5TFB2kUeOmxesk3Fkzy5y/Lh5wTp5NnBq166dNGzYUJYsWSKJiYlOfcCdO3dk0aJFUr9+fenQoYOrx0lkWIEBgRKRN0ItWCfyGNxvOSPMC+89YwkMFImIMC9YJ9049G3OmTNHlU/p3Lmzyt30wgsvyPfffy/Hjh2zuf3Ro0flu+++k+eee05t37VrV1XnDvshIiIiMnTmcK1e3YwZM+STTz5RgZFWjBSDxlGfLleuXHLz5k25du2aJCUlqdewa6QxQJ07BFGhoaFiFMwcbkzeyBwenxgvb6x5Q62/8+A7EhIU4vFjID+VGC+yx3zvSfV3RHjvGUd8vMgbydf2nXdEQnht9fr97XDgZA316JYtWyZ//vmn7NmzJ0Xm8OzZs6syLE2aNFFdfA888IC4CvsdP368at06efKk5M+fX9q0aSNjx46V4sWLu7zff/75R6VRQCmYBx98UFavXu30Phg4GRNLrpBfYckV42LJFbf9/nY4HYE1BEPWARFKsuAAcCA5dLo4CGpatGghW7ZsUUk3O3bsKFFRUTJ79mwVtOH5MmXKuLTvfv36qRY0IiIiImfoMmIMwVKxYsV0C5pg3LhxKjjCoPQjR47I/PnzZevWrfLBBx/IxYsXVbkXV3z55Zeybt06lfGciIiIyBk+OdQ+Pj5epk2bZin3goLBmsGDB6tutvXr18uOHTuc2i8GqA8dOlRatmwpPXr00P24iYiIyNh8MnDauHGj6vorW7as1KxZM83rmKUHS5cudWq/gwYNUuOmMMidiIiIyBCBE+rhQa1atWy+rj2PgenOJN9Ed9/rr78u5cqV0+lIiYiIyJ/4ZOCEGXRQokQJm69rz584ccKh/WHwOlIiVKxYUV577TUdj5SIiIj8iUuz6twtGlMnRSQ8PNzm69ogdOSNcsSbb76pgqy1a9dKiIu5LDALz3omHqYzEukBZVb2vbDPsk7kMSiz0tZ877HkisGgzMq+5GvLkivGD5z0tH37dpk6dao89dRT0qxZM5f3g3xSb731lq7HRgQos1K1cFVvHwb5I5RZyct7z5BQZqUqr63fdNVps+hiY2Ptdr0BspVnVCcPaQeQ2XzSpEmZOqYRI0aoAevacurUqUztj4iIiPykxenpp59WJVQaNGig/xGJSKlSpdTj6dOnbb6uPV+6dOl094Ptdu3aperldevWLcVrKA0DSGmgtUQhv5M9KBdjpJIx5DtQcuXdP99V6683eZ0lV8izJVf2m+89qfo6S64YreTKu8nX9vXXWXLF24HTV199JV9//bVUqVJFteg8+eSTki9fPt0OCiVbYOfOnTZf155HPidHnD9/Xi22IIBCTigib0lITJC31pu7gYfeP5SBE3mOKUFkX/IQhCpDRYT3nmEkJIhow0uGDmXg5O2uum+//VaVXNm/f7+88sorqm4cgifUsNNDo0aNVPmWY8eOqRaj1BYsWKAeO3TokO5+IiIiVKFhWwsGigNq1WnPEREREekeOPXs2VMFHiiWi0zcCHLmzp0rzZs3l8qVK6uyKJcuXRJXYebbwIED1fqAAQMsY5pg8uTJKn9T06ZNpXbt2pbnkWm8UqVKaiwSERERkc8NDkdm7wkTJqiB0mgFat26tSWYQq6l7t27y5o1a1zaN1II1K9fXzZt2iTly5eXxx9/XI2pGjJkiBQqVEhmzZqVYnsEaocPH5Zz585l5pSIiIiI3DurLjg4WLp06aKycx8/fly1EqHe3I8//iitWrVSmbo//PBDu7PkbAkLC1OtWiNHjlT5nBYvXqxyMUVGRqoxTmXKlNHj0ImIiIgcFmDScXDP77//LjNnzlRBDpJFIuBBeZTNmzdLUlKSlCxZUpYtWybVqlWTrA4JMNFFidQEuXPn1nXfEcOX67o/clzUhHYe/8yY+BjJOd6cgiN6RLTkCDEneCVyuzsxIj8kF1F/LFokmPeeYWCIS3JqH0FS6eTE0SSZ/v2d6Ran//77T3XXoTutZcuWqh4cWpiQdPLs2bNqwDhaoZ5//nlVSuWll17K7EcSERERZZ10BGik+uWXX1Tr0vLlyyUhIUHlOOrRo4cKkBo3bpxie4x3mj59uhqDtGXLFr2OncgQwoLD5K9n/7KsE3lMYJhI67/urpNxhIWJ/PXX3XXybuCEaf5ILokACq1L/fr1U0kxCxQokOH7tDQARGQWFBgkdYvX9fZhkD8KDBIpwHvPkIKCROry2vpM4IQuuM6dO6vWpYceesjh9w0bNkzleyIiIiLym8AJ6QdQxsRZFSpUUAsRpSy5MmXLFLU+qMEgZg4nz5ZcOWy+96TiIJZcMVrJlSnJ13bQIGYO15FLg8Nff/31NHmUbJkzZ4706dPHlY8g8quSK8NWD1ML1ok8WnJl1zDzgnUyVsmVYcPMC9bJu4ETAqINGzZkuN3GjRtVXTsiIiIiI9AlAaY9SIIZhAFqRERERAbgtsAJM+6Q4RvlUYiIiIj8anB4ixYtUvwbeZxSP6e5c+eOHDt2TM6fP89ZdEREROR/gdO6dess6wEBASoowmJPtmzZpH379jJp0qTMHyURERFRVgqcUDZF64JDgd2uXbvKxIkTbW4bEhIiBQsWVMETERERkd8FTqVLl7asjx49WmrWrJniOSJyDcqsrO1tzqjPkivkUSiz8mByNQeWXDEWlFnRKnWw5Ir3E2AicCIi/UquNIto5u3DIH8tuVKE954hYUZ7M17bLJeOgIiIiMjvWpwCAwPVcuDAAVUyxZncTBhIjll2RGQbsoV/vuNztd6vdj/JFsSxgeQhSQkiR833npTrJxLIe88wkC388+Rr268fZmx5+4j8K3AqVaqUCoC0wd4lS5ZU/yYifWrVDVw5UK1H3hfJwIk8JyleZLv53pMykQycjFarbmDytY2MZODk6cApKioq3X8TERER+QOOcSIiIiJyEAMnIiIiIj276k6ePCmZgTFSRERERH4ROEVERLg8GJyz6oiIiMivAqcHHniAs+iIiIjI7wU7W+CXiPQVGhwqy3oss6wTeUxgqEjTZXfXyThCQ0WWLbu7Tt4tuUJE+gkODJZ2Fdp5+zDIHwUGixTnvWdIwcEi7Xht3YGz6oiIiIj0bHH6448/1GO9evUkLCzM8m9HYYwUEdkvuTJ371y13ut/vZg5nDxbciXKfO9JRC9mDjdayZW5yde2Vy9mDtdRgMlkMmW0EerUYXD4wYMHVa067d+OSkxMFKO5ceOG5MmTR65fvy65c+fWdd8Rw5fruj9yXNQEzzdtx8THSM7xOdV69IhoyRGSw+PHQH7qTozID+Z7Tx6LFgnmvWcYMTEiOZOvbXS0SA5eW71+fzvU4vTUU0+pQAkfZP1vIiIiIn/iUOA0Z86cdP9NRERE5A84OJyIiIjIk+kI/vvvPzl79qxaL1asmBQpUkSP3RIREREZo8UJY8qnTp2qBosjWKpTp45asF6+fHmZMmWKJCUl6Xu0RERERFmtxSkuLk46dOgga9asUQFUvnz5pHTp0paCwMeOHZPBgwfLsmXL1BLKrKVERETkr4HTu+++K6tXr5Zq1arJxIkTpXXr1ile/+2332To0KHy+++/q23feustvY6XyHBQZuWHrj9Y1ok8BmVWGpvvPZZcMRg0WPyQfG3ZeOH5PE6plS1bVq5evSr//POPFChQwOY2ly5dUt14efPmlX///VeMhnmcjMkbeZyIiCjr/P52aYwTBoI/+OCDdoMmKFiwoLRo0ULOnTvnykcQERERGaOrrnjx4hIfH5/hdgkJCWqwOBHZdyfpjiw6uEitd67cWRX9JfKIpDsip833npTobC76S8Zw547IouRr27mzuegv6cKlFqdevXqpgeEnTpywuw1ewzY9e/Z0+eBu3bolo0aNUl1+qJGHIKxPnz5y5swZh/dx584dGTNmjLRr107KlCkjuXLlUvvCzL/+/funew5EnhB3J04eW/CYWrBO5DFJcSIbHjMvWCfjiIsTeewx84J18m7g9Oabb6puOBTvnTVrlsSgJk4yrM+ePVuaNm2quvMQ+Lji9u3b6jPGjh0r0dHR0rFjRylZsqTad82aNR0eN4X9YHA6ChMXLVpU2rRpowazo8Xsk08+kerVq8v27dtdOkYiIiLyLw613aGlJjWMKT99+rT07dtXLUhJABg0rkE9u0qVKqn0BM4aN26cbNmyRRo2bKhm6eVMLlY4efJkGTJkiGp5WrduXYb7QevShg0bpH79+hJs1VSJwsMIACdMmCDPP/88gyciIiLSZ1ZdYGDmKrM4mwgTrUGFCxdWI9537typWpis1ahRQ/bs2aOCndq1a7t8XOjGQ9cdWqWuXbtmKWLsCM6qMyZvzKqLiY+RnOPNfxhEj4iWHCGsYk4ecidG5AfzvSePRYsE894zDPQEJTc4SHS0SA5eW4/OqkPgk5nFWRs3blQnhLQHqYMm6Nq1q3pcunSpZAZaxIKCgtRjSEhIpvZFRERExueTRX53796tHmvVqmXzde15tDq5Cg1t7733nhqT1bx5c8mePbvL+yIiIiL/4JPzE1G2BUqUKGHzde15Z2fEvfbaa6ogMZrpEHRh7FXlypXliy++0OGoiYiIyOgyHTjdvHlTBSB4tDdcCrPvnIFZdBAeHm7z9RzJfbX4TGcsXLgwxUB1zKj79ttv5d5773WoPh8WDYIvIj2EBIXI7I6zLetEHhMYItJg9t11Mg4MP5mdfG05FMU3Aqd9+/bJyy+/rGa2ZTS+HDPYfMHRo0ct5WB27Nghb7zxhhpcPnPmTOndu3e67x0/fjxr7pFbZAvKJpH3RXr7MMgfBWYTKcN7z5CyZROJ5LX1mTFOqFHXuHFjVcQX6QK0Fpvu3btLvXr1LNP+H3nkEXnqqaec3r+WeiA2Ntbm61reKMyIcwXKwSCXExJ03nPPPfLCCy/IqVOn0n3PiBEj1IB1bcloeyIiIjIelwIn5FhCNxmSUf7555/SpEkT9fzcuXNl8+bNsn//fhVYHThwQOVdclapUqXUI/JE2aI9X7p0ackMTEfs0KGDylC+atWqdLcNDQ1V0xatFyK9Sq4sP7JcLVgn8hjcb2eWmxfee8YrubJ8uXnBOnk3cEJLEwZV2+veKleunPz8889y8eJFGTlypNP7R54mQA4nW7TnMUYps9D6BDhWIm9AmZX237VXC0uukEehzMr69uaFJVeMBWNy27c3Lyy54v3A6cKFC1KlShXLv7OhLzW5vIkmb9680qxZM1m2bJnT+2/UqJFqDcJA7l27dqV5fcGCBeoRrUWZtX79evWInFFEREREugdO+fPnTzHDDP+2lx4AQZazkIxy4MCBan3AgAEpauGh6w+pBFALzzpr+LRp01R5F4xFsrZ8+XLZtGlTms/A+CkMDkfghHFOqGFHREREpPusOgwGtw6S7rvvPjWzbv78+Zaivpi5hhl32nglZ6GO3OrVq1XQU758eTWOCp+5detWKVSokCoubA2fd/jwYTl37lyK57dt26ZmwxUvXlwdJ1qyzp8/r1qyrly5ov79ww8/WAakExEREena4tSqVSuVjkALntBlhrFCb7/9tppZhyK8devWVbPPHnvsMVc+QhXnXbt2rRojhXxOixcvVp8XGRmpxjjZKjxsS5cuXWTw4MFSrFgxFUQhSMIjBpajdergwYOWwe1EREREmS7ymxrGHiHbdufOnVX6AW3AOIIktOJoWrZsKUuWLFEz0oyGRX6NiUV+ya+wyK9xsciv235/u9RVh4HUSAhprUWLFqpFCOkJrl69KhUqVEgxBomIiIgoq9O1Vh1KoXCQNZFzUGZl2sPTLOtEHoMyK3XM9x5LrhgMyqxMS762LLnie4ETCueePXtWrWMsUZEiRfTYLZHflFwZUG+Atw+D/LXkSgXee4aENEEDeG19ZnA4YGjU1KlTVZccgqU6deqoBeuYBTdlyhRJSkrS92iJiIiIslqLE3I4YSYdar0hgMqXL5+l/MnJkyfV4HHMZEPySyxGHBxOpJfEpET58+Sfar1JqSYSFBjk7UMif5GUKHLRfO9JoSYivPeMIzFR5M/ka4uZ40G8tl5tcXr33XdVjqWqVavKypUr5fLlyypFABbkU/rll1+kWrVqaqYdtiUi+27fuS3Nv2quFqwTeUzSbZE1zc0L1sk4UMmjeXPzYlXVg7wUOH377beqpAryLLVu3dpmnie0RmG63zfffKPDYRIRERFl0cAJA8EffPBBKVCggN1tkBATKQpSZ/ImIiIi8qvACeVL4uPjM9wuISFBDRYnIiIi8tvAqVevXqorzlZRXw1ewzY9e/bMzPERERERZe3ACQV40Q33wAMPqGK7MUjtngzrs2fPlqZNm6ruPK3oLxEREZFfpCOwVVAXaQhOnz4tffv2VQtSEgDKrWgCAgKkUqVKKj0BERERkV8ETlFRURluY13cV5NeVx4R3c0c/v5D71vWiTwmIJvIfe/fXSdjZQ5///276+TZwIkZwIncB/XphjYa6u3DIH+E2ohVeO8ZEurTDeW19amSK0RERET+Rpciv0SUuZIrO8/tVOu1itZiyRXybMmVq+Z7T/LVYskVo5Vc2Zl8bWvVYskVX2lx2rNnjzz33HNSpUoVlSUcC9aff/559RoRZQxlVup9UU8tLLlCHoUyK7/WMy8suWIsKLNSr555YckV3wicpkyZInXq1JEvvvhCDh06JDdv3lQL1j///HP1GrYhIiIi8uvAadWqVfLKK69ISEiIevz7779VGoJr167Jrl27ZMiQIRIaGiqDBw9WSTCJiIiI/DZwmjx5sgQHB8tvv/0mkyZNkho1aqhuuty5c0v16tVl4sSJ6rXAwED54IMP9D9qIiIioqwSOP31118qM/j9999vd5uGDRtKs2bNZOvWrZk5PiIiIqKsHTjFxsZKoUKFMtwO22BbIiIiIr8NnEqWLCmbN2+WO3fu2N0Gr2EbbEtERETkt3mcOnbsqMYu9enTR6ZOnSp58+ZN8fqNGzdk0KBBcvLkSTVQnIjsQ5mV0U1HW9aJPAZlVqqZ7z2WXDEYlFkZnXxtWXJFVwEmVOt1EurS1a1bV9Wwy5kzp7Rp00YiIiIs9el++eUXFTyhOPC2bdssBYCNBOeHAfHXr19Xg+L1FDF8ua77I8dFTWjn7UMgIiIf/v3tUotT/vz55Y8//lCJLpcvXy4//vhjmm3atWsnn332mSGDJiIiIvJPLpdcKV68uCxdulSOHz8uGzZskLNnz6rnixUrJo0bN5Z7771Xz+MkMqwkU5IcvHhQrVcuVFkCA1hCkjzElCRy3XzvSZ7KIrz3jCMpSeRg8rWtXFkkkNfWq4FTrVq1pGzZsqqlCQESgyQi191KuCXVPqmm1qNHREuOkBzePiTyF4m3RFaY7z15LFokmPeeYdy6JVIt+dpGR4vk4LXVi0sh6OHDhyUbB5sRERGRn3EpcCpfvrxcvnxZ/6MhIiIiMlrg9Mwzz8j69etVQV8iIiIif+FS4PTiiy9KZGSkKrvy4YcfytGjRyU+Pl7/oyMiIiLK6oPDg4KC1CNSQL366qtqsScgICDdDONEREREhg6cUEYFARERERGRP3EpcELGcCLSB8qsvNrQ3GrLkivkUSizUjm5x4AlV4wFM9+13iDOgveNBJhEpI+QoBCZ2Gqitw+D/FFQiEhN3nuGFBIiMpHX1qcDp6tXr6pHFPxlNx4REREZUaZysC9ZskRatWqlCv0WLFhQLbly5VLP/fzzz/odJZGBoeRK1LUotWCdyGNwv0VHmRfee8YruYJhNViwTt4NnDCbrk+fPtK5c2dZvXq1xMbGqkrDWLCO57p06aJSFmBbV926dUtGjRolFSpUkLCwMFUHD5975swZh/dx7do1mTdvnvTo0UOVhgkJCVHBXf369WXKlCmSkJDg8vER6VVy5d4p96oF60QeLbmy5F7zgnUyVskVlEPDgnXybuCEgGPOnDlStGhR+eSTT1RwcuXKFbVcv35dPv30U/XaN998o7Z1xe3bt6VFixYyduxYiY6Olo4dO6rZfLNnz5aaNWvKv//+69B+Jk2aJL169ZL58+dLvnz5VEBXr1492b17t7z88svqMxDsEREREbklcPr8888lPDxc/vzzT3nuueckd+7cltfQmtOvXz/1Wvbs2dW2rhg3bpxs2bJFGjZsKEeOHFGBz9atW+WDDz6QixcvqpYnR+TIkUOGDRumZgLu3LlTvv/+e1mzZo3s3btXSpUqJRs2bFCfRUREROSWwOn48ePy4IMPqq4ve/AatsG2zkIW8mnTpqn16dOnqzFUmsGDB0v16tVVyZcdO3ZkuK8RI0bIe++9p4Kk1PX2JkyYoNa/++47p4+RiIiI/I9LgVOhQoXUWKGMZMuWTQ0Yd9bGjRtVl1/ZsmVVt1xqXbt2VY9Lly6VzKhRo4Z6PHv2bKb2Q0RERP7BpcAJg8J///13SwoCWzDeCdt06tTJ6f1j/BHUqlXL5uva83v27JHM0MZJ3XPPPZnaDxEREfkHlwInjAkqU6aMGliN4Ci1tWvXSsuWLVWL0bvvvuv0/k+ePKkeS5QoYfN17fkTJ05IZmgD1zHwnIiIiMgtCTARaKCrDmOMECDlz59fSpcubQl6Ll++rNYbNGiQJihBckwMzk4PZtEBBqDbG/ANN2/eFFdh5h/SJiBh5/DhwzPcPi4uTi2aGzduuPzZRNaCA4Olf53+lnXKeiKGL5esKCQgQd4s2k6tjxu5SuJNWa80R9QE8/FTKsHBIv37310n3bj0ba5bt86yjjxNCJS0YMna5s2b0zznC1nFMeNv0KBB6lhmzZql8kNlZPz48fLWW2955PjIv4QGh8r0dtO9fRjkhxAojTr7grcPg9whNBSzq7x9FIbkUuDkykw5Z2iz6OzlV4qJibGkPnDWvn37VCsYZu5NnTpVjddyBGbnYUafdYsT8koRERGR/3ApcNK65dxFSx1w+vRpm69rzzt7HAj4UA4Gg9rHjBkjL774osPvDQ0NVQuR3tBqeyn2klovGF7QJ1plyV+YJH+QedjBlUTk4+O9Zxio2nHJ/HNFMLudP1d045Mdn1qaACSstEV7HvmcHHXu3Dk1HguP6KYbPXq0TkdLlDmxCbFSeFJhtR49IlpyhJjH8BG5W/aAONlZtZdar7x3gdwyhXn7kEgv6LEpbP65Ihg3nDw2mLxc5NddGjVqpOreHTt2THbt2pXm9QULFqjHDh06OLQ/tDC1bt1a7e/pp5+WDz/8UPdjJiIiIuPzycAJM/YGDhyo1gcMGGAZ0wSTJ09W+ZuaNm0qtWvXtjyPTOOVKlVSY5GsYZxUu3btVImVxx57TGbOnMmuECIiIjJOVx28+eabKl3Apk2bVHmUJk2aqLxNqFeHzOWYDWft0qVLcvjwYdUVZ+2NN95Qs/uCgoIkODhYnnnmGZufh6LFRERERFkycAoLC1OJNJEGYN68ebJ48WKVLyoyMlLGjh1rNzlmalp288TERLUfexg4ERERUZbsqtNkz55d3n77bTl69KhKPonWpNmzZ9sMmjBLDrOTUgdA+Deez2ghIiIiytKBExEREZEv8dmuOiJ/gTIrvWv0tqwTeUqiBMmCKw9a1slAUGalt/nnCkuu6IvfJpEPlFyZ04lj7Mg7JVdePf2Ktw+D3AEJmzl21y3YVUdERETkILY4EXkZJicgeziEZwtnnjHyIJPKHg63TCgpxXvPMDDpSav3Gh7Okis6YosTkZchaMo5PqdatACKyBMQNB38X1e1aAEUGQSCppw5zYsWQJEuGDgREREROYiBExEREZGDGDgREREROYiBExEREZGDGDgREREROYiBExEREZGDmMeJyMuCAoOka5WulnUiT0mSQFl+rZFlnQwkKEika9e766QbBk5EXhYWHCY/dvvR24dBfijOFCIDTo7w9mGQO4SFifzInyvuwD8xiIiIiBzEwImIiIjIQQyciLwsJj5GAt4KUAvWiTwle8BtiareXi1YJwOJiTHXp8OCddINAyciIiIiBzFwIiIiInIQAyciIiIiBzFwIiIiInIQAyciIiIiBzFwIiIiInIQM4cTeRnKrLQt39ayTuQpKLPy+406lnUyEJRZaWv+ucKSK/pi4ETkAyVXlvdc7u3DID8tudInaoy3D4PcVXJlOX+uuAP/xCAiIiJyEAMnIiIiIgcxcCLyMpRZyfFuDrWw5Ap5EsqsHKj2qFpYcsVgUGYlRw7zwpIruuIYJyIfEJsQ6+1DID8VHhjn7UMgd4nlzxV3YIsTERERkYMYOBERERE5iIETERERkYMYOBERERE5iIETERERkYM4q47IywIDAqVp6aaWdSJPSZIA2RJdzbJOBhIYKNK06d110g0DJyIvy54tu6yLXOftwyA/FGcKle7/TvD2YZA7ZM8uso4/V9yBYSgRERGRgxg4ERERERkhcLp165aMGjVKKlSoIGFhYVKsWDHp06ePnDlzxqn9rF+/Xt566y1p166dFCpUSAICAiQiIsJtx03kDJRZKTSxkFpYcoU8CWVWdlTpqRaWXDEYlFkpVMi8sOSKf4xxun37trRo0UK2bNkiRYsWlY4dO0pUVJTMnj1bli1bpp4vU6aMQ/saNGiQ7N692+3HTOSqS7GXvH0I5KcKBN/w9iGQu1zizxW/anEaN26cCo4aNmwoR44ckfnz58vWrVvlgw8+kIsXL6qWJ0e1atVK7e/XX3+V/fv3u/W4iYiIyLh8ssUpPj5epk2bptanT58uOXPmtLw2ePBg+eqrr1T3244dO6R27doZ7u/999+3rJ8/f95NR02UsYjhy9M8lyS3RbKb1yuP+kUCJczzB+YHoia08/YhEJEB+GSL08aNG+X69etStmxZqVmzZprXu3btqh6XLl3qhaMjIiIif+WTgZM2HqlWrVo2X9ee37Nnj0ePi4iIiPybTwZOJ0+eVI8lSpSw+br2/IkTJzx6XEREROTffHKMU3R0tHoMDw+3+XqOHDnU482bNz12THFxcWrR3LjBmSiklwAJSSpvWSfyFJRZ2R1rvvdYcsVgUGalTp2762TswMkXjR8/XuWCItJboIRK0bgPvX0Y5KclVzoe5b1n2JIr27Z5+ygMySfDUG0WXWxsrM3XY5KTeeXKlctjxzRixAg1YF1bTp065bHPJiIiIt/gky1OpUqVUo+nT5+2+br2fOnSpT12TKGhoWohIiIi/+WTLU41atRQjzt37rT5uvZ89erVPXpcRO6APE6nQ/uoReV0IvKQsIDbsqFSH7VgnQwEPTYoLYbFTu8NGajFqVGjRpInTx45duyY7Nq1S+67774Ury9YsEA9dujQwUtHSKSvxMAL3j4E8kMYDl4ixHzvcWi4wZhMmHp+d52M3eIUEhIiAwcOVOsDBgywjGmCyZMnq/xNTZs2TZE1HJnGK1WqpMYiEREREflNixO8+eabsnr1atm0aZOUL19emjRpovI2oV5doUKFZNasWSm2v3Tpkhw+fFjOnTuXZl9ffPGFWiAhIUE9YrsGDRpYtpkxY4bdhJtEREREPh04hYWFydq1a1UagHnz5snixYslf/78EhkZKWPHjrWbHNPeYHIEXKnr4Vk/x7xMRERElCW76jTZs2eXt99+W44ePaqST6KVaPbs2TaDpjFjxojJZJI5c+bYfS29pVmzZh46KyIiIsqqfDpwIiIiIvIlPttVR+RPsiWZc5cReRLmWh25bb73OO/KYAICRKpUubtOumHgRORlgRImxeJmePswyA/dNoVJqyO89wwJtV737/f2URgSu+qIiIiIHMTAiYiIiMhBDJyIvAxlVs6G9lcLS66QJ6HMym8V+quFJVcMBmVWqlY1Lyy5oiuOcSLyAQmBJ719COSHMGS4Qpj53uPwYYNBmZUDB+6uk27Y4kRERETkIAZORERERA5i4ERERETkIAZORERERA5i4ERERETkIM6qI/IBQUmFvX0I5Icw1+p0vPne47wrg0GZldKl766Tbhg4EflAyZUScbO8fRjkpyVXGh/ivWfYkitRUd4+CkNiVx0RERGRgxg4ERERETmIgRORlyVJnJwLfUUtWCfylNCAOPm53CtqwToZyK1bInXrmhesk244xonI60wSH/iPZZ3IUwLFJDXC/7Gsk4EkJYls3353nXTDFiciIiIiB7HFiYiIyI6I4cslK8oef1sOJq9XHvmL3AoJk6wmakI78UVscSIiIiJyEAMnIiIiIgcxcCIiIiJyEMc4EfmAQFNubx8C+anLd3jvGdXl7Ly27sDAicgHSq6UvD3P24dBfuiWKUxqH+C9Z0QYDF77JV5bd2BXHREREZGDGDgREREROYiBE5GXoczK+ZDhamHJFfIklFn5vsxwtbDkirGEJsTJ9/OGqwXrpB+OcSLyOpPEBe2zrBN5CsqsNMhpvvdYcsVYAk0maXAq+dqaeG31xBYnIiIiIgcxcCIiIiJyEAMnIiIiIgcxcCIiIiJyEAMnIiIiIgdxVh2RDwgwhXr7EMhPxSbx3jOq2Gy8tu7AwInIB0qulLq90NuHQX5acqXKPt57Ri25UmUwr607sKuOiIiIyEEMnIiIiIiMEDjdunVLRo0aJRUqVJCwsDApVqyY9OnTR86cOeP0vq5evSqDBg2S0qVLS2hoqHp8+eWX5dq1a245diJHmSReLoSMUQvWiTwlNCBeZkWMUQvWyThC78TLrB/HqAXr5AdjnG7fvi0tWrSQLVu2SNGiRaVjx44SFRUls2fPlmXLlqnny5Qp49C+Ll26JA0bNpSjR4+q93Tq1En2798vU6ZMkZUrV8rmzZslf/78bj8nIltMkiS3grZb1gO8fUDkNwIlSVrk3m5ZJ+MITEqSFv8mX9skXlu/aHEaN26cCo4Q8Bw5ckTmz58vW7dulQ8++EAuXryoWp4chZYlBE1dunSRw4cPq33t27dPXnzxRbXvwYMHu/VciIiIyBh8MnCKj4+XadOmqfXp06dLzpw5La8hyKlevbqsX79eduzYkeG+zp07J999952EhITIjBkzJDj4biPbxIkTpVChQvLtt9/KhQsX3HQ2REREZBQ+GTht3LhRrl+/LmXLlpWaNWumeb1r167qcenSpRnu65dffpGkpCRp0qSJFClSJMVrGOvUoUMHSUxMlBUrVuh4BkRERGREPhk47d69Wz3WqlXL5uva83v27PHovoiIiMi/+WTgdPLkSfVYokQJm69rz584ccKj+yIiIiL/5pOz6qKjo9VjeHi4zddz5MihHm/evOmxfcXFxalFg65EuHHjhugtKS5W932SY9xxPTO6tklyW7SpdObXOQMmq13brPr/bGLAbbmRfOiJcbGSZMp69543/p/NChLjb4v2zfDa2t6vyWQSwwROvmj8+PHy1ltvpXm+ZMmSXjkeco88H3n388/IU949AAPz9rX1VXksa1nz3uN1deDazuC1tQUNJnny3P0/IEsHTtosuthY25F+TEyMesyVK5fH9jVixIgUaQsw4PzKlStSoEABCQgISDeyRXB16tQpyZ07txidP50vz9W4/Ol8ea7G5U/ne8OJc0VLE4ImJNV2hU8GTqVKlVKPp0+ftvm69jyyf3tqX5iBh8Va3rx5xVG4kEa/cf31fHmuxuVP58tzNS5/Ot/cDp6rKy1NPj04vEaNGupx586dNl/Xnkc+J0/ui4iIiPybTwZOjRo1UtHgsWPHZNeuXWleX7BggXpEDqaMtGnTRgIDA+XPP/9Mk+QSg72RCyooKEjatm2r4xkQERGREflk4IQs3wMHDlTrAwYMsIxDgsmTJ6ucS02bNpXatWtbnkem8UqVKqmxSNZQ565Hjx4qG3n//v3lzp07lteGDRumyrc88cQTUrhwYbecC7r3Ro8enaabz6j86Xx5rsblT+fLczUufzrfUA+ea4DJ1fl4Hijy26xZM1WfDsEPMn8j1xL+jTIpqYv8jhkzRs166927t8yZMydNkd8GDRqoFixkI69Tp44q8ot6deXLl1f7YpFfIiIiypItThAWFiZr166VkSNHqhxMixcvVoFTZGSkGpdkHTRlpGDBgvLXX3+por5oeVq0aJHKw/TSSy+p5xk0ERERUZZucSIiIiLyNT7b4kRERETkaxg46Wzjxo1qhh66/5B8s169evL11187vR+M00JiTXtL9+7dxd1u3bolo0aNkgoVKqiuUyQL69Onj5w5c8bpfV29elUGDRqk8mVh8B4eX375Zbl27Zr4Cr3ONyIiIt1rd+jQIfGmHTt2yIQJE6RLly6qVqN2XK7y5Wur57n6+nVFkl8MaXjmmWekYsWK6h5GSSmkZHn77bct5aeMcG31Pldfv7baxCjcxxiXi1nn2vV46qmnZO/evYa5tnqfqzuuLbvqdLRw4UJ5/PHHVVbxBx54QI2tWrNmjboRhwwZIpMmTXIqcHr66afVD4L77rsvzev169eXF154Qdw5OL958+Zq4Lw2OD8qKkqNCbM1OD89GJzfsGFDOXr0qHqPNjgfC4KUzZs3e32cmZ7ni/9RMR4PExXsle/BZ3hLp06d5Oeff07zvCs/Cnz92up5rr5+Xb/44gvp27evWq9cubJUq1ZNZVPetGmTypKMWcfr1693eAaxL19bvc/V168t4PcJZpgj52Dx4sXVc7gWR44ckWzZsslPP/0k7du3z/LXVu9zdcu1ReBEmXf58mVT7ty58dPYtHDhQsvz58+fN5UrV049v3btWof3N3v2bPWe0aNHm7zhjTfeUJ/fsGFD082bNy3Pf/DBB+r5pk2bOryvXr16qfd06dLFlJCQYHn+xRdfVM/37t3b5G16nm/p0qXVe3zVhAkTTCNHjjQtWbLEdO7cOVNoaKjLx+vr11bPc/X16zpnzhxTv379TAcOHEjx/NmzZ001a9ZUx96jRw9DXFu9z9XXry1s2LDBdOvWrTTPT58+XR17kSJFUlynrHpt9T5Xd1xb375TspD33ntPXZyOHTumee2nn35Sr7Vv3z5LBE5xcXGmPHnyqM/fuXNnmterV6+uXtu+fXuG+8IPssDAQFNISIgKIq3dvn3bVKhQIVNQUJDpv//+M3mLnuebVX4IW3M1mMgK1zY1IwdO6dm0aZM6dpw/7ncjXltXzzWrX1soW7asOv7du3cb+to6e67uurYc46ST5cuXq8euXbumea1du3aqD3716tWqSygrjNNCugbkvKpZs2aa17VzRNb1jPzyyy+q6xJdX0WKFEnxGvqtkf09MTFRVqxYIUY4X3+SFa4tpSw9hWoJly9fNvS1dfZcjQDdV1ryaCNfW2fP1V18sshvVrR79271WKtWrTSv4QKjD3779u2qj9aZungY3Dp06FDVf3/PPfdIixYtVNZ0b52L9fPI4K7HvmbNmuXQvrLC+VqbOHGiSrqKH0ZVq1aVzp07q/FSRpEVrq07ZMXr+u+//1p+6TgydiUrX1tnzzWrX9tvvvlGDh8+rAZSYzHytf3GyXN117Vl4KQDBDVosQDM3LEFzyNwwiA1ZwKnZcuWqUWDGSMInObPn5/mrwW9nDx50nLMtmjP41w8uS93cdcxoqSPtVdeeUU+/vhjNVPPCLLCtXWHrHhdp0yZYqnd6UhJiqx8bZ0916x2bREAYKA0Bk8fPHhQrWMG8Hfffafqrhrp2k7M5Lm669qyq04H1lNfkeXcFkyVBcz4cARG+aOMzN9//62CsvPnz8uSJUsss0UwowDNqe48Hz3ORc99uYvex/jII4+oWR/4wYNp0yjtM3jwYNV18Oyzz9qc5ZUVZYVrq6esel3R5fLll1+qFpixY8ca+tq6cq5Z7dr++uuv8tVXX6li9wgkME0fgYR17VajXNtfM3mubru2uo6YysI6depkqlixolPL1q1b1XvPnDmjBp9hsTfSX5vFMHfu3EwdJ2Z8VahQQe1r3rx5Jnfo27ev2j9mmtnyzz//qNfLly+f4b5atmyptp05c6bN11etWqVex3beouf5pufzzz9X+8G940tcHTCdFa6tnoPDs9p1hYMHD5ry5cunju+jjz4y9LV19Vyz6rW9evWq6Y8//jA99NBD6hjHjRtn2Gt71cVzdde1ZYtTsuPHj6u+U2cWRK+ARJca7bnU0NQIuXLlytRx4rNQY0+Lxt1BOx89zkXPfbmLp44RyfqQVwb3DnJEZXVZ4dp6gq9eVyRuRXcVEh3iL2wkOzTqtc3MuWbFawt58+ZVA7zRyoYWGNR13bZtm+GubWbO1V3XloFTsl27dqmkeM4szZo1U+/NnTu3ym4Kp0+ftrl/7Xk0NWaWNiju3Llz4g6lSpXS7Vz03Je7eOoYAwMD1cw9d147T8oK19YTfPG6XrlyRVq1aqW6J5BI15nku1nt2mb2XLPatU0N3ZJIvIzfSY7M/M1K1zaz5+qua8vASecpsDt37kzzWkJCgupXRUoCZGXNLPxVZd0X7clzsX7ekUHueu7LXTx5jO6+dp6UFa6tp/jSdcUYlocfflgOHDigylbMnDnT6TIzWeXa6nGuWenappdpGy5evGiYa6vHubrt2ma6o5DckgAzPd26dVP7Gzt2rMndCSH//vtv3RJgpk6o5ivJ1vQ83/Ts27fPFBAQYAoPD3c4MV9WSYDpq9fWE2OcfOm64ntv0aKFOsfWrVu7fDxZ4drqda5Z5dqmB5m+8T1MnDjRENdWr3N117Vl4OTmkiu4+dIruaINND99+nSK5999913TxYsXUzwXHx9vGjNmjNpX9uzZ07zHHSVI7r//flN0dHSGJUg+/vhjdR7Dhw+3OzD+0UcfTTF4/qWXXvKJ9P56nu/y5ctNa9asSbN/ZLmtXLmy2hfO25dkFExk9Wurx7lmhet6584dU+fOndWxNGnSxBQTE5Phe7LqtdXzXLPCtUUJkpUrV5oSExPT/E6YOnWqCoTwO+HkyZNZ/tpu0PFc3XVtGTjpaMGCBeqiIopt3ry5qWvXrqa8efOqizN48GCb79Fm4x0/fjzN8/gh36hRI1P37t1Nbdu2NRUrVkw9HxYWliI4cwfUCapfv776vKJFi5oee+wxy7/xF8mxY8dSbI/SMPb+h0MAqKXJx+Pjjz9uqlatmmWmGoJOb9PrfLXnkeb/kUceUdeuXr16puDgYPV8s2bNTLGxsSZvWrZsmTo3bcH9imOzfg7bGOHa6nWuWeG6YiaZ9vMEQQXOwdZi/QdZVr22ep5rVri2WgmuggULqta1nj17mlq1aqV+Vmm/E+bPn5/iPVn12s7W8VzddW0ZOLkhWm7Tpo0KmNAEWKdOHVWQ0h57gdOoUaPUdNBSpUqp6Bo3C1qunnvuOdOhQ4c8cCYmdUOhQCr+x0Kz7j333GOKjIw0nTp1Ks226f1PCvgfEQUkS5YsqfaFR0T6mGbqK/Q4X9TJ6tOnj+l///ufqUCBAup/0Pz586v/QTH9F38pe5v2gym9BdsY4drqda5Z4bpqx57RYv2zJqteWz3PNStc23///df0+uuvqz+kEUBky5bNlCNHDlPVqlXV9UHKlNSy6rX9V8dzdde1DcB/Mj3CioiIiMgPcFYdERERkYMYOBERERE5iIETERERkYMYOBERERE5iIETERERkYMYOBERERE5iIETERERkYMYOBERERE5iIETEekOFer1rlLvTZGRkep85syZo8v+7ty5IxUrVpR69eqlu926deukWbNm6W5z7tw5yZ49u/Tv31+XYyOi9DFwIiLysM8++0yOHDkiY8aMyfS+ihYtKv369ZOZM2eqfRKRezFwIiLyoLi4OHn77belevXq0rZt2zSvR0dHy7Bhw6Rs2bLSqlUrWb9+veTPn19t/9xzz8nRo0fTvAfbJyUlyciRIz10FkT+i4ETEZEHLViwQC5cuCBPPfVUmtcQ/LRu3VomTpwo165dk/r160uRIkWkadOmcvPmTfn8889l165dad5XvHhxad68uSxatEj+++8/D50JkX9i4EREXnfq1CnVmlK6dGkJDQ2VwoULS5cuXWTbtm02t9+3b5888cQTUqZMGQkLC5NChQrJfffdJy+//LIa82Nt06ZN0qlTJ8u+77nnHjW2aPjw4ap1x9O++OILNV6qe/fuaV5buXKlOt4aNWqolqWxY8dKpUqVVEB0/PhxNeapcuXKNvfbs2dPSUhI0G0cFhHZxsCJiLxq7969UqtWLdWagkHOCJjKly+vgoX7779ffvzxxxTb79ixQ+rWrStz586VXLlySceOHaVBgwYqaJgyZYocPnzYsu3SpUulSZMmsmTJEjUWCPuuWbOmXLlyRd577z25dOmSR8/1xo0b8ueff0q5cuVUK1Fqe/bsUY9PP/205MuXL83raHmqWrWqzX1rg8iXL1+u+3ET0V3BVutERB5lMpmkV69eKoDBOJ0JEyZYZuMtXLhQHnvsMenTp480btxYBT4wdepUuX37tkyaNEmGDBmSYn+HDh2SPHnyWP6NbdD9he6xRx99NMW2aM0qUKCAeBJakxITE1XgZ0vOnDnV4+nTp53eN1rfChYsKH/99Zf6ftASR0T6Y4sTEXkNup7Q4lSqVCkZN25cihQGCHTQxYbutFmzZlmev3jxonp86KGH0uwP3VpagJXRtghe0GLlSVqLElIR2ILB4sHBwSo4fOutt+Sff/5xav/YLwafHzx4UJfjJaK0GDgRkdeg2wrQspQtW7Y0rz/55JMptoPatWurxwEDBqjACzmR7NG2xX7QwoTWJ2/CoHCw1Q0HmEn3zTffSEhIiEpVgDQDW7duVbPrPvzwQzVgPD2YfWcdMBKR/hg4EZHXnD17Vj1GRETYfF17/syZM5bnhg4dqsbzbNy4Uc0kQxCCwALjm65fv57i/e+++64aaI2xThgQjq6sRx55RA3QRneWp2nHl15LFwaNR0VFyfTp06Vly5aqBWnVqlUyePBgqVChgurusyd37tzqMaMAi4hcx8CJiHyWrezjCA5+//131QqFcVFVqlRR/8aMOnRVWXdvlSxZUrZv3y6//vqrvPjii+rfCKL69u2r8iJdvnzZo+ejjb9CaoH0YOwVMoG//vrr8sADD6jElkhfgJYkBFYYJ5VeYJY3b143HD0RAQMnIvKaYsWKqccTJ07YfB0tL5B6BhoCKgwYx8w4dGWh5apHjx4qh9Ebb7yRYluMGUKLFMYN7d69W+2zRYsWKsDC+z0JaRYAs/qcgVmGSDOA2YdI3YB0DLZcvXpVPSI9AxG5BwMnIvIapAoApByw1Yry7bffptguvYBEK19iL6jQIJ/Ta6+95tC2ekO3IVinTHAUgkW0mAG672zBrELkqrKX64mIMo+BExF5DcYq/e9//1OtQKNGjVLpCTTI4/TTTz+pKfpISaD59NNPVTLI1FasWKEeteACMKD6/PnzDm3rCchLFRQUZDex53fffae6FW1Ba9nq1atVritbuZyOHTumuh4xloupCIjch3mciMhtkJjSnmeffVYtSGSJQd4YyI1gCRnAT548qQZ/o5vtyy+/TJFiAIHTCy+8oMY2oWUF26ClBYEFAgYEYBpM6X/11VdVSw+6uxCYYTuMGcIMNLzmDGTyxufbgmPE8acHg8LReobZgMjVVKJEiRSvoyUKx4wEmehOjI+PV0EiEncuW7ZMJflErqscOXKk2Tf2Ce3atXPqnIjISSYiIp3hR0tGy+jRoy3bnzhxwtS3b19TyZIlTdmyZTMVLFjQ1KlTJ9PWrVvT7HvJkiWmPn36mKpWrWrKmzevKTw83FShQgXTs88+azp06FCKbb/++mtTz549TRUrVjTlypVLLVWqVDENHjzYdPr0aYfPp3fv3hmeT+nSpR3a19y5c9X277//fprXLl++bJoyZYqpdevWpjJlyphCQ0PVtjlz5jQ1atTING/ePLv7bdGihfruzp8/7/B5EZHzAvAfZ4MtIiJyDcYnYZwVxmVpCTHtQSsSxm5prUn2oPUK++zatavMnz9f5yMmImsc40RE5EEYvI3uRGRMR/ebHiZOnCiBgYHy9ttv67I/IrKPgRMRkYchIziSWWI8U2adO3dOFUhGbip7pVyISD/sqiMiIiJyEFuciIiIiBzEwImIiIjIQQyciIiIiBzEwImIiIjIQQyciIiIiBzEwImIiIjIQQyciIiIiBzEwImIiIjIQQyciIiIiBzEwImIiIhIHPN/MajQGUEwiMUAAAAASUVORK5CYII=",
"text/plain": [
" ┌───────┐┌────────┐ ┌───────────┐\n",
" state_0: ┤0 ├┤0 ├──────┤0 ├\n",
" │ ││ │ │ │\n",
" state_1: ┤1 ├┤1 ├──────┤1 ├\n",
" │ P(X) ││ │ │ │\n",
" state_2: ┤2 ├┤2 ├──────┤2 ├\n",
" │ ││ │ │ │\n",
" state_3: ┤3 ├┤3 ├──────┤3 ├\n",
" └───────┘│ adder │┌────┐│ adder_dg │\n",
"objective: ─────────┤ ├┤2 ├┤ ├\n",
" │ ││ ││ │\n",
" sum_0: ─────────┤4 ├┤0 F ├┤4 ├\n",
" │ ││ ││ │\n",
" sum_1: ─────────┤5 ├┤1 ├┤5 ├\n",
" │ │└────┘│ │\n",
" carry: ─────────┤6 ├──────┤6 ├\n",
" └────────┘ └───────────┘"
],
"text/plain": [
" ┌───────┐┌────────┐ ┌───────────┐\n",
" state_0: ┤0 ├┤0 ├──────┤0 ├\n",
" │ ││ │ │ │\n",
" state_1: ┤1 ├┤1 ├──────┤1 ├\n",
" │ P(X) ││ │ │ │\n",
" state_2: ┤2 ├┤2 ├──────┤2 ├\n",
" │ ││ │ │ │\n",
" state_3: ┤3 ├┤3 ├──────┤3 ├\n",
" └───────┘│ adder │┌────┐│ adder_dg │\n",
"objective: ─────────┤ ├┤2 ├┤ ├\n",
" │ ││ ││ │\n",
" sum_0: ─────────┤4 ├┤0 F ├┤4 ├\n",
" │ ││ ││ │\n",
" sum_1: ─────────┤5 ├┤1 ├┤5 ├\n",
" │ │└────┘│ │\n",
" carry: ─────────┤6 ├──────┤6 ├\n",
" └────────┘ └───────────┘"
]
},
"execution_count": 15,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# define the registers for convenience and readability\n",
"qr_state = QuantumRegister(u.num_qubits, \"state\")\n",
"qr_sum = QuantumRegister(agg.num_sum_qubits, \"sum\")\n",
"qr_carry = QuantumRegister(agg.num_carry_qubits, \"carry\")\n",
"qr_obj = QuantumRegister(1, \"objective\")\n",
"\n",
"# define the circuit\n",
"state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name=\"A\")\n",
"\n",
"# load the random variable\n",
"state_preparation.append(u.to_gate(), qr_state)\n",
"\n",
"# aggregate\n",
"state_preparation.append(agg.to_gate(), qr_state[:] + qr_sum[:] + qr_carry[:])\n",
"\n",
"# linear objective function\n",
"state_preparation.append(objective.to_gate(), qr_sum[:] + qr_obj[:])\n",
"\n",
"# uncompute aggregation\n",
"state_preparation.append(agg.to_gate().inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])\n",
"\n",
"# draw the circuit\n",
"state_preparation.draw()"
]
},
{
"cell_type": "code",
"execution_count": 16,
"id": "ede38a99-cd81-41fc-a86f-c2c79778ffea",
"metadata": {},
"outputs": [],
"source": [
"state_preparation_measure = state_preparation.measure_all(inplace=False)\n",
"sampler = Sampler(backend = qr_backend)\n",
"job = sampler.run([state_preparation_measure])\n",
"\n",
"from qiskit.result import QuasiDistribution\n",
"result = job.result()\n",
"pub_result = result[0]\n",
"counts = pub_result.data.meas.get_counts()\n",
"shots = 0\n",
"for key, val in counts.items():\n",
" shots += val\n",
"quasi_dist = QuasiDistribution({outcome: freq / shots for outcome, freq in counts.items()})\n",
"\n",
"binary_probabilities = quasi_dist.binary_probabilities()"
]
},
{
"cell_type": "code",
"execution_count": 17,
"id": "25185c53-6b89-4619-8009-1400cf164590",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact Expected Loss: 0.6592\n",
"Exact Operator Value: 0.3867\n",
"Mapped Operator value: 0.6346\n"
]
}
],
"source": [
"# evaluate the result\n",
"value = 0\n",
"for i, prob in binary_probabilities.items():\n",
" if prob > 1e-6 and i[-(len(qr_state) + 1) :][0] == \"1\":\n",
" value += prob\n",
"\n",
"print(\"Exact Expected Loss: %.4f\" % expected_loss)\n",
"print(\"Exact Operator Value: %.4f\" % value)\n",
"print(\"Mapped Operator value: %.4f\" % objective.post_processing(value))"
]
},
{
"cell_type": "code",
"execution_count": 18,
"id": "d1067c98-f6f5-4f44-8ced-a40754b6e9fa",
"metadata": {},
"outputs": [],
"source": [
"# Example output\n",
"\n",
"# Exact Expected Loss: 0.6562\n",
"# Exact Operator Value: 0.0000\n",
"# Mapped Operator value: -2.3197"
]
},
{
"cell_type": "code",
"execution_count": 19,
"id": "b45b6004-e91a-424a-9bc3-a4f2ffec2b2f",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact value: \t0.6592\n",
"Estimated value:\t0.6517\n",
"Confidence interval: \t[0.5830, 0.7204]\n"
]
}
],
"source": [
"# set target precision and confidence level\n",
"epsilon = 0.01\n",
"alpha = 0.05\n",
"\n",
"problem = EstimationProblem(\n",
" state_preparation=state_preparation,\n",
" objective_qubits=[len(qr_state)],\n",
" post_processing=objective.post_processing,\n",
")\n",
"# construct amplitude estimation\n",
"ae = IterativeAmplitudeEstimation(\n",
" epsilon_target=epsilon, alpha=alpha, sampler=Sampler(backend = qr_backend, options={\"shots\": 100, \"seed\": 75})\n",
")\n",
"result = ae.estimate(problem)\n",
"\n",
"# print results\n",
"conf_int = np.array(result.confidence_interval_processed)\n",
"print(\"Exact value: \\t%.4f\" % expected_loss)\n",
"print(\"Estimated value:\\t%.4f\" % result.estimation_processed)\n",
"print(\"Confidence interval: \\t[%.4f, %.4f]\" % tuple(conf_int))"
]
},
{
"cell_type": "code",
"execution_count": 20,
"id": "97a187a6-3ba1-4e09-9acc-9813af63db4d",
"metadata": {},
"outputs": [],
"source": [
"# Example output:\n",
"\n",
"# Exact value: \t0.6562\n",
"# Estimated value:\t0.6592\n",
"# Confidence interval: \t[0.5982, 0.7202]"
]
},
{
"cell_type": "code",
"execution_count": 21,
"id": "81933d68-67e1-4ae0-a5d1-df12c92e8331",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" ┌──────────────┐\n",
"state_0: ┤0 ├\n",
" │ │\n",
"state_1: ┤1 ├\n",
" │ circuit-148 │\n",
"compare: ┤2 ├\n",
" │ │\n",
" a17: ┤3 ├\n",
" └──────────────┘"
],
"text/plain": [
" ┌──────────────┐\n",
"state_0: ┤0 ├\n",
" │ │\n",
"state_1: ┤1 ├\n",
" │ circuit-148 │\n",
"compare: ┤2 ├\n",
" │ │\n",
" a17: ┤3 ├\n",
" └──────────────┘"
]
},
"execution_count": 21,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# set x value to estimate the CDF\n",
"x_eval = 2\n",
"\n",
"comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)\n",
"comparator.draw()"
]
},
{
"cell_type": "code",
"execution_count": 22,
"id": "a6358b35-7bdf-479d-acc2-868d84fb3c0c",
"metadata": {},
"outputs": [],
"source": [
"def get_cdf_circuit(x_eval):\n",
" # define the registers for convenience and readability\n",
" qr_state = QuantumRegister(u.num_qubits, \"state\")\n",
" qr_sum = QuantumRegister(agg.num_sum_qubits, \"sum\")\n",
" qr_carry = QuantumRegister(agg.num_carry_qubits, \"carry\")\n",
" qr_obj = QuantumRegister(1, \"objective\")\n",
" qr_compare = QuantumRegister(1, \"compare\")\n",
"\n",
" # define the circuit\n",
" state_preparation = QuantumCircuit(qr_state, qr_obj, qr_sum, qr_carry, name=\"A\")\n",
"\n",
" # load the random variable\n",
" state_preparation.append(u, qr_state)\n",
"\n",
" # aggregate\n",
" state_preparation.append(agg, qr_state[:] + qr_sum[:] + qr_carry[:])\n",
"\n",
" # comparator objective function\n",
" comparator = IntegerComparator(agg.num_sum_qubits, x_eval + 1, geq=False)\n",
" state_preparation.append(comparator, qr_sum[:] + qr_obj[:] + qr_carry[:])\n",
"\n",
" # uncompute aggregation\n",
" state_preparation.append(agg.inverse(), qr_state[:] + qr_sum[:] + qr_carry[:])\n",
"\n",
" return state_preparation\n",
"\n",
"\n",
"state_preparation = get_cdf_circuit(x_eval)"
]
},
{
"cell_type": "code",
"execution_count": 23,
"id": "3e901ba0-5e81-41f3-99d8-b93ae8c0d084",
"metadata": {},
"outputs": [
{
"data": {
"text/html": [
" ┌───────┐┌────────┐ ┌───────────┐\n",
" state_0: ┤0 ├┤0 ├────────┤0 ├\n",
" │ ││ │ │ │\n",
" state_1: ┤1 ├┤1 ├────────┤1 ├\n",
" │ P(X) ││ │ │ │\n",
" state_2: ┤2 ├┤2 ├────────┤2 ├\n",
" │ ││ │ │ │\n",
" state_3: ┤3 ├┤3 ├────────┤3 ├\n",
" └───────┘│ adder │┌──────┐│ adder_dg │\n",
"objective: ─────────┤ ├┤2 ├┤ ├\n",
" │ ││ ││ │\n",
" sum_0: ─────────┤4 ├┤0 ├┤4 ├\n",
" │ ││ cmp ││ │\n",
" sum_1: ─────────┤5 ├┤1 ├┤5 ├\n",
" │ ││ ││ │\n",
" carry: ─────────┤6 ├┤3 ├┤6 ├\n",
" └────────┘└──────┘└───────────┘"
],
"text/plain": [
" ┌───────┐┌────────┐ ┌───────────┐\n",
" state_0: ┤0 ├┤0 ├────────┤0 ├\n",
" │ ││ │ │ │\n",
" state_1: ┤1 ├┤1 ├────────┤1 ├\n",
" │ P(X) ││ │ │ │\n",
" state_2: ┤2 ├┤2 ├────────┤2 ├\n",
" │ ││ │ │ │\n",
" state_3: ┤3 ├┤3 ├────────┤3 ├\n",
" └───────┘│ adder │┌──────┐│ adder_dg │\n",
"objective: ─────────┤ ├┤2 ├┤ ├\n",
" │ ││ ││ │\n",
" sum_0: ─────────┤4 ├┤0 ├┤4 ├\n",
" │ ││ cmp ││ │\n",
" sum_1: ─────────┤5 ├┤1 ├┤5 ├\n",
" │ ││ ││ │\n",
" carry: ─────────┤6 ├┤3 ├┤6 ├\n",
" └────────┘└──────┘└───────────┘"
]
},
"execution_count": 23,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"state_preparation.draw()"
]
},
{
"cell_type": "code",
"execution_count": 24,
"id": "8deaaa33-4088-4bac-a556-8d7761d8316b",
"metadata": {},
"outputs": [],
"source": [
"state_preparation_measure = state_preparation.measure_all(inplace=False)\n",
"sampler = Sampler(backend = qr_backend)\n",
"job = sampler.run([state_preparation_measure])\n",
"\n",
"from qiskit.result import QuasiDistribution\n",
"result = job.result()\n",
"pub_result = result[0]\n",
"counts = pub_result.data.meas.get_counts()\n",
"shots = 0\n",
"for key, val in counts.items():\n",
" shots += val\n",
"quasi_dist = QuasiDistribution({outcome: freq / shots for outcome, freq in counts.items()})\n",
"\n",
"binary_probabilities = quasi_dist.binary_probabilities()"
]
},
{
"cell_type": "code",
"execution_count": 25,
"id": "9c0a54e9-a23b-4b45-a9e3-3761ac02cf18",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Operator CDF(2) = 0.9668\n",
"Exact CDF(2) = 0.9512\n"
]
}
],
"source": [
"# evaluate the result\n",
"var_prob = 0\n",
"for i, prob in binary_probabilities.items():\n",
" if prob > 1e-6 and i[-(len(qr_state) + 1) :][0] == \"1\":\n",
" var_prob += prob\n",
"\n",
"print(\"Operator CDF(%s)\" % x_eval + \" = %.4f\" % var_prob)\n",
"print(\"Exact CDF(%s)\" % x_eval + \" = %.4f\" % cdf[x_eval])"
]
},
{
"cell_type": "code",
"execution_count": 26,
"id": "80bb90b4-413a-4985-b41d-0d1e95dff5a3",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Exact value: \t0.9512\n",
"Estimated value:\t0.9783\n",
"Confidence interval: \t[0.9775, 0.9790]\n"
]
}
],
"source": [
"# set target precision and confidence level\n",
"epsilon = 0.01\n",
"alpha = 0.05\n",
"\n",
"problem = EstimationProblem(state_preparation=state_preparation, objective_qubits=[len(qr_state)])\n",
"# construct amplitude estimation\n",
"ae_cdf = IterativeAmplitudeEstimation(\n",
" epsilon_target=epsilon, alpha=alpha, sampler=Sampler(backend = qr_backend, options={\"shots\": 100, \"seed\": 75})\n",
")\n",
"result_cdf = ae_cdf.estimate(problem)\n",
"\n",
"# print results\n",
"conf_int = np.array(result_cdf.confidence_interval)\n",
"print(\"Exact value: \\t%.4f\" % cdf[x_eval])\n",
"print(\"Estimated value:\\t%.4f\" % result_cdf.estimation)\n",
"print(\"Confidence interval: \\t[%.4f, %.4f]\" % tuple(conf_int))"
]
},
{
"cell_type": "code",
"execution_count": 27,
"id": "a52476a6-280b-40a8-adf8-fde3aa188498",
"metadata": {},
"outputs": [],
"source": [
"# Example output:\n",
"\n",
"# Exact value: \t0.9600\n",
"# Estimated value:\t0.9581\n",
"# Confidence interval: \t[0.9570, 0.9592]"
]
},
{
"cell_type": "code",
"execution_count": 28,
"id": "fd99ddcc-dd14-479e-ad17-db16d500f9fe",
"metadata": {},
"outputs": [],
"source": [
"def run_ae_for_cdf(x_eval, epsilon=0.01, alpha=0.05):\n",
"\n",
" # construct amplitude estimation\n",
" state_preparation = get_cdf_circuit(x_eval)\n",
" problem = EstimationProblem(\n",
" state_preparation=state_preparation, objective_qubits=[len(qr_state)]\n",
" )\n",
" ae_var = IterativeAmplitudeEstimation(\n",
" epsilon_target=epsilon, alpha=alpha, sampler=Sampler(backend = qr_backend, options={\"shots\": 100, \"seed\": 75})\n",
" )\n",
" result_var = ae_var.estimate(problem)\n",
"\n",
" return result_var.estimation"
]
},
{
"cell_type": "code",
"execution_count": 29,
"id": "fa19ad5e-3487-4991-9c3d-616cd30970f6",
"metadata": {},
"outputs": [],
"source": [
"def bisection_search(\n",
" objective, target_value, low_level, high_level, low_value=None, high_value=None\n",
"):\n",
" \"\"\"\n",
" Determines the smallest level such that the objective value is still larger than the target\n",
" :param objective: objective function\n",
" :param target: target value\n",
" :param low_level: lowest level to be considered\n",
" :param high_level: highest level to be considered\n",
" :param low_value: value of lowest level (will be evaluated if set to None)\n",
" :param high_value: value of highest level (will be evaluated if set to None)\n",
" :return: dictionary with level, value, num_eval\n",
" \"\"\"\n",
"\n",
" # check whether low and high values are given and evaluated them otherwise\n",
" print(\"--------------------------------------------------------------------\")\n",
" print(\"start bisection search for target value %.3f\" % target_value)\n",
" print(\"--------------------------------------------------------------------\")\n",
" num_eval = 0\n",
" if low_value is None:\n",
" low_value = objective(low_level)\n",
" num_eval += 1\n",
" if high_value is None:\n",
" high_value = objective(high_level)\n",
" num_eval += 1\n",
"\n",
" # check if low_value already satisfies the condition\n",
" if low_value > target_value:\n",
" return {\n",
" \"level\": low_level,\n",
" \"value\": low_value,\n",
" \"num_eval\": num_eval,\n",
" \"comment\": \"returned low value\",\n",
" }\n",
" elif low_value == target_value:\n",
" return {\"level\": low_level, \"value\": low_value, \"num_eval\": num_eval, \"comment\": \"success\"}\n",
"\n",
" # check if high_value is above target\n",
" if high_value < target_value:\n",
" return {\n",
" \"level\": high_level,\n",
" \"value\": high_value,\n",
" \"num_eval\": num_eval,\n",
" \"comment\": \"returned low value\",\n",
" }\n",
" elif high_value == target_value:\n",
" return {\n",
" \"level\": high_level,\n",
" \"value\": high_value,\n",
" \"num_eval\": num_eval,\n",
" \"comment\": \"success\",\n",
" }\n",
"\n",
" # perform bisection search until\n",
" print(\"low_level low_value level value high_level high_value\")\n",
" print(\"--------------------------------------------------------------------\")\n",
" while high_level - low_level > 1:\n",
"\n",
" level = int(np.round((high_level + low_level) / 2.0))\n",
" num_eval += 1\n",
" value = objective(level)\n",
"\n",
" print(\n",
" \"%2d %.3f %2d %.3f %2d %.3f\"\n",
" % (low_level, low_value, level, value, high_level, high_value)\n",
" )\n",
"\n",
" if value >= target_value:\n",
" high_level = level\n",
" high_value = value\n",
" else:\n",
" low_level = level\n",
" low_value = value\n",
"\n",
" # return high value after bisection search\n",
" print(\"--------------------------------------------------------------------\")\n",
" print(\"finished bisection search\")\n",
" print(\"--------------------------------------------------------------------\")\n",
" return {\"level\": high_level, \"value\": high_value, \"num_eval\": num_eval, \"comment\": \"success\"}"
]
},
{
"cell_type": "code",
"execution_count": 30,
"id": "c1933cc0-c5ab-4df7-aea5-24e7e123dd86",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"--------------------------------------------------------------------\n",
"start bisection search for target value 0.950\n",
"--------------------------------------------------------------------\n",
"low_level low_value level value high_level high_value\n",
"--------------------------------------------------------------------\n",
"-1 0.000 1 0.751 3 1.000\n",
" 1 0.751 2 0.958 3 1.000\n",
"--------------------------------------------------------------------\n",
"finished bisection search\n",
"--------------------------------------------------------------------\n"
]
}
],
"source": [
"# run bisection search to determine VaR\n",
"objective = lambda x: run_ae_for_cdf(x)\n",
"bisection_result = bisection_search(\n",
" objective, 1 - alpha, min(losses) - 1, max(losses), low_value=0, high_value=1\n",
")\n",
"var = bisection_result[\"level\"]"
]
},
{
"cell_type": "code",
"execution_count": 31,
"id": "2721d65b-624d-4610-88ac-634db3a8c8c1",
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"Estimated Value at Risk: 2\n",
"Exact Value at Risk: 2\n",
"Estimated Probability: 0.958\n",
"Exact Probability: 0.951\n"
]
}
],
"source": [
"print(\"Estimated Value at Risk: %2d\" % var)\n",
"print(\"Exact Value at Risk: %2d\" % exact_var)\n",
"print(\"Estimated Probability: %.3f\" % bisection_result[\"value\"])\n",
"print(\"Exact Probability: %.3f\" % cdf[exact_var])"
]
},
{
"cell_type": "code",
"execution_count": 32,
"id": "1a334390-8df7-494c-baa1-9ab2905321b7",
"metadata": {},
"outputs": [],
"source": [
"# Example output:\n",
"\n",
"# Estimated Value at Risk: 2\n",
"# Exact Value at Risk: 2\n",
"# Estimated Probability: 0.960\n",
"# Exact Probability: 0.960"
]
},
{
"cell_type": "code",
"execution_count": 33,
"id": "8729fb66-bdfc-42a1-9024-c9d2eb25c2d5",
"metadata": {},
"outputs": [
{
"name": "stderr",
"output_type": "stream",
"text": [
"C:\\Users\\vkasi\\AppData\\Local\\Temp\\ipykernel_39140\\2011604251.py:9: DeprecationWarning: The class ``qiskit.circuit.library.arithmetic.linear_amplitude_function.LinearAmplitudeFunction`` is deprecated as of Qiskit 2.2. It will be removed in Qiskit 3.0. Use the class qiskit.circuit.library.LinearAmplitudeFunctionGate instead.\n",
" cvar_objective = LinearAmplitudeFunction(\n"
]
},
{
"data": {
"text/html": [
" ┌────┐\n",
"q78_0: ┤0 ├\n",
" │ │\n",
"q78_1: ┤1 ├\n",
" │ │\n",
" q79: ┤2 F ├\n",
" │ │\n",
"a80_0: ┤3 ├\n",
" │ │\n",
"a80_1: ┤4 ├\n",
" └────┘"
],
"text/plain": [
" ┌────┐\n",
"q78_0: ┤0 ├\n",
" │ │\n",
"q78_1: ┤1 ├\n",
" │ │\n",
" q79: ┤2 F ├\n",
" │ │\n",
"a80_0: ┤3 ├\n",
" │ │\n",
"a80_1: ┤4 ├\n",
" └────┘"
]
},
"execution_count": 33,
"metadata": {},
"output_type": "execute_result"
}
],
"source": [
"# define linear objective\n",
"breakpoints = [0, var]\n",
"slopes = [0, 1]\n",
"offsets = [0, 0] # subtract VaR and add it later to the estimate\n",
"f_min = 0\n",
"f_max = 3 - var\n",
"c_approx = 0.25\n",
"\n",
"cvar_objective = LinearAmplitudeFunction(\n",
" agg.num_sum_qubits,\n",
" slopes,\n",
" offsets,\n",
" domain=(0, 2**agg.num_sum_qubits - 1),\n",
" image=(f_min, f_max),\n",
" rescaling_factor=c_approx,\n",
" breakpoints=breakpoints,\n",
")\n",
"\n",
"cvar_objective.draw()"
]
},
{
"cell_type": "code",
"execution_count": 34,
"id": "1a0d1e86-3190-4a2b-8bdf-e7594b169afb",
"metadata": {},
"outputs": [
{
"data": {
"text/plain": [
"| Software | Version |
|---|---|
qiskit | 2.2.1 |
qiskit_algorithms | 0.4.0 |
qiskit_finance | 0.4.1 |
qiskit_aer | 0.17.1 |
qiskit_ibm_runtime | 0.40.1 |
| System information | |
| Python version | 3.12.9 |
| OS | Windows |
| Thu Oct 16 11:28:31 2025 Mountain Daylight Time | |
| Software | Version |
|---|---|
QuantumRingsLib | 0.11.0 |
quantumrings-toolkit-qiskit | 0.1.10 |
© Copyright IBM 2017, 2025.
This code is licensed under the Apache License, Version 2.0. You may
obtain a copy of this license in the LICENSE.txt file in the root directory
of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.
Any modifications or derivative works of this code must retain this
copyright notice, and modified files need to carry a notice indicating
that they have been altered from the originals.
Modifications Copyright Quantum Rings Inc, 2025
Modified from the originals
Added support for Quantum Rings QrEstimatorV2 class.