docs/tutorials/03_european_call_option_pricing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pricing European Call Options "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "<br>\n",
    "Suppose a <a href=\"http://www.theoptionsguide.com/call-option.aspx\">European call option</a> with strike price $K$ and an underlying asset whose spot price at maturity $S_T$ follows a given random distribution.\n",
    "The corresponding payoff function is defined as:\n",
    "\n",
    "$$\\max\\{S_T - K, 0\\}$$\n",
    "\n",
    "In the following, a quantum algorithm based on amplitude estimation is used to estimate the expected payoff, i.e., the fair price before discounting, for the option:\n",
    "\n",
    "$$\\mathbb{E}\\left[ \\max\\{S_T - K, 0\\} \\right]$$\n",
    "\n",
    "as well as the corresponding $\\Delta$, i.e., the derivative of the option price with respect to the spot price, defined as:\n",
    "\n",
    "$$\n",
    "\\Delta = \\mathbb{P}\\left[S_T \\geq K\\right]\n",
    "$$\n",
    "\n",
    "The approximation of the objective function and a general introduction to option pricing and risk analysis on quantum computers are given in the following papers:\n",
    "\n",
    "- [Quantum Risk Analysis. Woerner, Egger. 2018.](https://www.nature.com/articles/s41534-019-0130-6)\n",
    "- [Option Pricing using Quantum Computers. Stamatopoulos et al. 2019.](https://quantum-journal.org/papers/q-2020-07-06-291/)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:09.974086Z",
     "start_time": "2020-07-13T23:35:09.967567Z"
    }
   },
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from qiskit import QuantumCircuit\n",
    "from qiskit_algorithms import IterativeAmplitudeEstimation, EstimationProblem\n",
    "from qiskit.circuit.library import LinearAmplitudeFunction\n",
    "from qiskit_aer.primitives import Sampler\n",
    "from qiskit_finance.circuit.library import LogNormalDistribution"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uncertainty Model\n",
    "\n",
    "We construct a circuit to load a log-normal random distribution into a quantum state.\n",
    "The distribution is truncated to a given interval $[\\text{low}, \\text{high}]$ and discretized using $2^n$ grid points, where $n$ denotes the number of qubits used.\n",
    "The unitary operator corresponding to the circuit implements the following: \n",
    "\n",
    "$$\\big|0\\rangle_{n} \\mapsto \\big|\\psi\\rangle_{n} = \\sum_{i=0}^{2^n-1} \\sqrt{p_i}\\big|i\\rangle_{n},$$\n",
    "\n",
    "where $p_i$ denote the probabilities corresponding to the truncated and discretized distribution and where $i$ is mapped to the right interval using the affine map:\n",
    "\n",
    "$$ \\{0, \\ldots, 2^n-1\\} \\ni i \\mapsto \\frac{\\text{high} - \\text{low}}{2^n - 1} * i + \\text{low} \\in [\\text{low}, \\text{high}].$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:10.294975Z",
     "start_time": "2020-07-13T23:35:10.283925Z"
    }
   },
   "outputs": [],
   "source": [
    "# number of qubits to represent the uncertainty\n",
    "num_uncertainty_qubits = 3\n",
    "\n",
    "# parameters for considered random distribution\n",
    "S = 2.0  # initial spot price\n",
    "vol = 0.4  # volatility of 40%\n",
    "r = 0.05  # annual interest rate of 4%\n",
    "T = 40 / 365  # 40 days to maturity\n",
    "\n",
    "# resulting parameters for log-normal distribution\n",
    "mu = (r - 0.5 * vol**2) * T + np.log(S)\n",
    "sigma = vol * np.sqrt(T)\n",
    "mean = np.exp(mu + sigma**2 / 2)\n",
    "variance = (np.exp(sigma**2) - 1) * np.exp(2 * mu + sigma**2)\n",
    "stddev = np.sqrt(variance)\n",
    "\n",
    "# lowest and highest value considered for the spot price; in between, an equidistant discretization is considered.\n",
    "low = np.maximum(0, mean - 3 * stddev)\n",
    "high = mean + 3 * stddev\n",
    "\n",
    "# construct A operator for QAE for the payoff function by\n",
    "# composing the uncertainty model and the objective\n",
    "uncertainty_model = LogNormalDistribution(\n",
    "    num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEyCAYAAADOV2anAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5FklEQVR4nO3deZgcVbnH8e9PEAiEJbKERSSACiJRMEGMoCSArHqDIERFr0E0ogJeRRARIYALoCxeuIoRNSJqVEQUBDEsCSJ7FAlLlC3EgIJAIGaDhLz3j3MGOpXume6Znqom+X2ep5+ZqjpV9XZPT79ddTZFBGZmZv3tFVUHYGZmKwcnHDMzK4UTjpmZlcIJx8zMSuGEY2ZmpXDCMTOzUlSecCRtJ+laSQskPSbpVEmr9LDPGyX9Ppd/TtIsSRdK2qRQbqKkqPPYtn+flZmZFa1a5cklDQKuAe4FRgNbA2eREuGJ3ey6LvAwcBHwGLAlcDIwTNJOEbGkpuwM4LDC/jObiW+DDTaIIUOGNFO0X8yfP5+11lqrsvM30qlxgWPrjU6NCxxbb1Qd17Rp056MiA3rboyIyh7AF4E5wDo1644DFtSua/JY7wICeEvNuonAHb2Nb9iwYVGl66+/vtLzN9KpcUU4tt7o1LgiHFtvVB1Xd5+5Vd9S2xe4OiLm1qybBAwAdmvxWE/ln6u1IzAzM2uvqhPOtqRbXi+KiFmkK5we61kkvULSapK2AU4HbgduKxTbTtLcXNdzo6RWE5mZmbWBosKx1CQtBo6NiHML62cDF0XECT3s/3tg77w4DdgvIp6o2f4Z4HlSHdGGwDHAMGDXiCgmpq59xgHjAAYPHjxs0qRJvXhm7TFv3jwGDhxY2fkb6dS4wLH1RqfGBY6tN6qOa9SoUdMiYnjdjY3utZXxABYD/1Nn/Wzga03s/zpgZ+BDpCulacAa3ZRfk9TY4LJm4nMdTn2dGleEY+uNTo0rwrH1RtVx0cF1OHNILc6KBuVt3YqI+yPi1oi4mHSlsyPwwW7KLwCuBN7Su3DNzKy3qk44MyjU1UjanHQlMqPuHg1ExCPA08BWPRXNDzMzK1HVCecqYG9Ja9esGwMsBKa2cqDccGB90i2zRmUGAPuTbr2ZmVmJKu34CVwAHA1cKukM0tXJeODsqGkqLekBYGpEHJ6XvwksAW4FngHeQOq/8yCpWTWS1gWuAC4GHgA2AD4LbAoc3P9PzczMalWacCJijqQ9gPOBy0nJ4xxS0qm1KlA73M0dwFGk1mRrALOAXwFfj4j5ucxzwL9JIxZsBCwCbgZ2i4g7+uHpmJlZN6q+wiEi7gV276HMkMLyJPKVTDf7LAIO7Gt8ZgBDjv9dn49xzNAljO3jcWaevn+f4zCrStV1OGZmtpJwwjEzs1I44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZlcIJx8zMSuGEY2ZmpXDCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSVJ5wJG0n6VpJCyQ9JulUSav0sM8bJf0+l39O0ixJF0rapE7Z0ZKmS1ok6V5JY/rv2ZiZWSOrVnlySYOAa4B7gdHA1sBZpER4Yje7rgs8DFwEPAZsCZwMDJO0U0QsycffFfgV8G3gaGA/4GeS5kTEH/rlSZmZWV2VJhzgCGAAcGBEzAUmS1oHGC/pzLxuORFxE3BTzaopkmYDfwDeBPw5r/8ycENEHJ2Xr5f0RuCkXNbMzEpS9S21fYGrC4llEikJ7dbisZ7KP1cDkLQ6MAr4RaHcJGCEpHVbD9fMzHqr6oSzLTCjdkVEzAIW5G3dkvQKSatJ2gY4HbgduC1v3hp4ZfH4wH2k5/36voVuZmatUERUd3JpMXBsRJxbWD8buCgiTuhh/98De+fFacB+EfFE3rYLcCOwY0TcWbPPa4H7gb3r1eNIGgeMAxg8ePCwSZMm9e7JtcG8efMYOHBgZedvpFPjgv6Lbfqjz/b5GIMHwOML+3aMoZu1/8J8Zfx7tkOnxlZ1XKNGjZoWEcPrbau6DqevjgJeBbyO1MjgKkm7RMSi3h4wIiYAEwCGDx8eI0eObEecvTJlyhSqPH8jnRoX9F9sY4//XZ+PcczQJZw1vW//cjMPHdnnOIpWxr9nO3RqbJ0aF1SfcOaQWpwVDcrbuhUR9+dfb5X0R1LLtQ8CP6jZv3j8QTXnNjOzklRdhzODQl2NpM2BNVm+7qVbEfEI8DSwVV71ILC4ePy8vBT4ey/iNTOzXqo64VwF7C1p7Zp1Y4CFwNRWDpQbDqxPusohIp4DrgcOLhQdA9wcEX2/KW9mZk2r+pbaBaQOmZdKOoN0dTIeOLu2qbSkB4CpEXF4Xv4msAS4FXgGeANwHOmqpraW/zRSH51zgctIHT/3A/bpx+dkZmZ1VJpwImKOpD2A84HLScnjHFLSqbUqUDvczR2kBgPjgDWAWaQRBb4eEfNrjn+jpPcBXwE+Sa7j8SgDtqIY0qbGDH1tFDHz9P37HIet+Kq+wiEi7gV276HMkMLyJJa9kulu38tIVzdmZlahqutwzMxsJeGEY2ZmpXDCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwULY8WLWko8FZgY9LUAE+TZs+8KSI8bbOZmdXVVMKRtBVpPplDgcGkKZqfAZ4D1iNNCb1U0lTgQuDnEbG0H+I1M7OXqR5vqUm6ELgH2AE4FdgRWCMiNoyIV0fEQGAj4D3AdOBM4D5Ju/Zb1GZm9rLTzBXOQmDbiHikUYGIeBK4CrhK0ueAg4HN2hOimZmtCHq8womIo7pLNnXKL42In0fEz5spL2k7SddKWiDpMUmnSlqlh312kvRDSQ/k/f4m6WRJaxTKjZcUdR77NPt8zMysPfo0xbSk7YHdAAFTI2J6i/sPAq4B7gVGA1sDZ5ES4Ynd7Domlz0DuB94E3Ba/nlQoeyzQDHB3NdKnGZm1ne9TjiSPgl8FbgWWAv4hqRjIuLbLRzmCGAAcGBEzAUmS1oHGC/pzLyuntPzbbwuUyQtAr4raYvCFdmSiLilhZjMzKwfNNNoYM0Gm74AjIiIgyNiP+BI4Estnn9f4OpCYplESkK7NdqpkGy6/CX/3LTFGMzMrATNdPz8u6RD66wXqXl0l940g94WmFG7IiJmAQvytlaMyDE8WFi/nqQnJS2W9BdJB/YiTjMz6yNFRPcFpHcC5wLPA0dHxG15/adJzaSvJfXD2QM4LiLOa/rk0mLg2Ig4t7B+NnBRRJzQ5HE2Bu4CroyIsTXrP0Rqsv0XYG3gE8B+wEERcWmDY40DxgEMHjx42KRJk5p9Om03b948Bg4cWNn5G+nUuKD/Ypv+6LN9PsbgAfD4wr4dY+hm6y6z3KlxtcvK+F7rq6rjGjVq1LSIGF5vW48JB0CSgMNJFfOTgS9ExD8lvZmXbn3dEBF3thJYOxKOpNVIDQ9eDQzrbrSD/DxuAgZExA49HXv48OFxxx139FSs30yZMoWRI0dWdv5GOjUu6L/Yhhz/uz4f45ihSzhrep/a6TDz9P2XWe7UuNplZXyv9VXVcUlqmHCaGkstkguBbYDHgbslfQmYERH/mx939iK2OUC9r0aD8rZu5QRyEfBGYL+ehtaJlF0vBd7UU9NrMzNrr5YG74yIuRFxLLAzaTy1GZLe14fzz6BQVyNpc9Ituhl191jWuaTm1KMjopnyAJEfZmZWoqZaqUn6iqRbc6X7BGBRRIwm1XWcLGlqvr3WqquAvSWtXbNuDGl0g6k9xPVFUsu4D0XEjc2cLF8RHQT8NSJe6EW8ZmbWS83cuP0+sB2pz80CUpKZLGm7iJicE82n8rrLImJcC+e/ADgauFTSGcBWwHjg7Nqm0pIeIHUsPTwvfxD4GjAReFTS22qO+WBE/DuXmwr8inS1tBbwcdLV2QEtxGhmZm3QTMLZFzg4IiYDSPoT8BSpp/8DeVTo8yX9FDi5lZNHxBxJewDnA5eTRqA+h5R0inHW1rnslX+OzY9ah5ESEcADwP8Am5CaTP8Z2D8irmolTjMz67tmEs4M4MOSpgGLSE2L5wOzawtFxNPAZ1oNICLuBXbvocyQwvJYlk809fY7vNV4zMysfzSTcD5CumJ4klTZPpN0xbOo/8IyM7MVTY8JJyL+BoyQtBawmmf1NDOz3mi6t1dEzCfdSjMzM2tZM82iP9xqJ0lJr5X0jt6HZWZmK5pmOn5+DnhQ0mnd9bWRtL6kQyVdDtxJahlmZmYGNFeHs6OkMcBRwJckzSNNYPYk8BywHrAl8BrScDQXA0dExKP9FbSZmb38NFWHk6eL/rmkrYE9gbcAG5M6Uz4O3AD8CZgSEYv7KVYzM3sZa2mI2Ih4kOXnmzEzM+tRS4N3mpmZ9ZYTjpmZlcIJx8zMSuGEY2ZmpWgp4Uh6jyQnKTMza1mryeMyYLakMyS9oR/iMTOzFVSrCWdr4HvAIcDdkm6W9HFJ67Q/NDMzW5G0lHAiYmZEnBwRWwLvIk1wdg7wT0k/ljSqP4I0M7OXv17Xx0TEdRHxYeD1wDTgUOAaSQ9J+qykljqVmpnZiq3XCUfSbpImAn8Dtgf+jzT18yXAKcBF7QjQzMxWDK22UttC0kmSHgSuAzYHxgGbRMRREXFtRBxHmiV0dJPH3E7StZIWSHpM0qk9TYcgaSdJP5T0QN7vb5JOlrRGnbK7SLpV0iJJD0s6upXnbGZm7dHqba+HgMdIU07/ICIeblDuHuC2ng4maRBwDXAvKUFtDZxFSoQndrPrmFz2DOB+4E3AafnnQTXHfy1wNXAF8EXgrcDZkhZExIU9xWdmZu3TasJ5N3B1RCztrlBE/B1opgHBEcAA4MCImAtMzi3exks6M6+r5/SIeLJmeYqkRcB3JW0REY/k9ceSEuSHImIJcJ2k1wAnS/p+REQTMZqZWRu0WoezE2laguVI2kTSSS0eb19SAqtNLJNISWi3RjsVkk2Xv+SfmxaOf2lONrXHfzWp3snMzErSasI5mfRhXc+meXsrtgVm1K6IiFnAgrytFSOApeTpEyStRapjmlEod1/Nuc3MrCRq5a6SpKXAzhFxe51to4HvR8QGLRxvMXBsRJxbWD8buCgiTmjyOBsDdwFXRsTYvG4zYDbw3oi4rKbsqsBi4BMRMaHOscaRGkIwePDgYZMmTWr26bTdvHnzGDhwYGXnb6RT44L+i236o8/2+RiDB8DjC/t2jKGbrbvMcqfG1S4r43utr6qOa9SoUdMiYni9bT3W4Uj6CKnVGUAA35FUrFtZAxgK/KEvgfaGpNWAXwDzgM/29Xg5CU0AGD58eIwcObKvh+y1KVOmUOX5G+nUuKD/Yht7/O/6fIxjhi7hrOl9654289CRyyx3alztsjK+1/qqU+OC5hoNLACeyr8LeBZ4ulDmeeAq4Nstnn8OUO+r0aC8rVuSROrv80Zgl4io3eeZ/LN4/EE15zYzs5L0mHAi4pfALwEk/RA4tZvm0K2aQaEuRdLmwJosX/dSz7mk5tTviohiXdB8Sf8oHr9muZnjm5lZm7Q6ltphbUw2kK6K9pa0ds26McBCYGp3O0r6InAkqcnzjd0c/72FjqRjgH8Ad/c6ajMza1nVc9tcADwHXCppz1xhPx44u7apdB5R4Ps1yx8Evka6nfaopLfVPDasOf43SK3qfixplKTjgE+QrtLcB8fMrETNNBq4DRgbEfdKup3UcKChiHhrsyePiDmS9gDOBy4n1bucQ0o6xThrr1L2yj/H5ketw0gjIRARD0jaBzibdLXzL+AYjzJgZla+ZhoN3EO6xdX1e1uvDCLiXmD3HsoMKSyPZflE02jfG0lD2piZWYWaaTRwWM3vY/s1GjMzW2FVXYdjZmYriWbqcHqst6nVSh2OmZmtPJqtw3GLLjMz65Nm6nDGlhCHmZmt4FyHY2Zmpai0H46Zma08Ku+HY2ZmKwf3wzEzs1K0PAlGnn9mLKn3/ibAP4FbgR9FxPNtjc7MzFYYLTUakPQG4H7g/4DtgRfyz/8DHpC0XdsjNDOzFUKrVzgTSBOwvSMiZnWtlPQa4ArS6M/vbF94Zma2omg14QwHPlCbbAAiYpakk4Gfti0yW+kMadN0yX2ddnnm6fv3OQ4zW16r/XBmAms02LYGMKvBNjMzW8m1mnCOB74iaefalZLeBpwGfKFdgZmZ2YqlN4N3rgPcJOkJ4Algo/x4CjgBuKz9YZqZ2ctdbwbvvKefYjEzsxVY5YN35qbU5wEjSFNMXwicEhEvdLPPasBXgbeRGjKsERGqU24i8JE6h3hDRMzoc/BmZta0ljt+tpOkQcA1wL3AaGBr4CxS3dKJ3ey6JvAx4DbgJrqfonoGcFhh3czeRWxmZr1VacIBjgAGAAdGxFxgsqR1gPGSzszrlhMRz0h6VUSEpCPpPuHMj4hb2h+6mZm1ouXpCSSNkXSNpFmSnig+WjzcvsDVhcQyiZSEdutux4jwIKJmZi8jrQ5t80HgR8ADwKuB35JGGHgFMBc4v8Xzb0u65fWi3Kl0Qd7WDttJmivpOUk3Suo2kZmZWf9QKxcKkv4CXAKcDiwGhkfEnyWtDUwGLomIb7ZwvMXAsRFxbmH9bOCiiDihiWMcCZzXoNHAZ4DnSXVEGwLHAMOAXSPitgbHGweMAxg8ePCwSZMmNft02m7evHkMHDiwsvM30l9xTX/02T4fY/AAeHxhz+W6M3SzdZdb16mxdWpc7dKp/wPQubFVHdeoUaOmRcTwettarcN5HfCniHhB0gukPjlExH8knQGcAzSdcPpbRHyrdlnSlaRm3ScABzTYZwJpzDiGDx8eI0eO7N8guzFlyhSqPH8j/RVXX4ekgTS0zVnT+1Y1OfPQkcut69TYOjWudunU/wHo3Ng6NS5ovQ5nLrB6/v1R4A012wSs3+Lx5gD1vhoNytvaKiIWAFcCb2n3sc3MrHutfq25HXgTcDWp/uYkSUtIt61OAlptDTaDQl2NpM1JzZ77q59M4FlLzcxK12rC+TqwRf79pPz7d0hXSrcDn2jxeFcBx0paOyL+k9eNIU1pPbXFY/VI0gBgf2Bau49tZmbdaynh5P4st+TfnwFGS1odWL1Rn5keXAAcDVya64C2AsYDZ9ceT9IDwNSIOLxm3b7AWsAOefl9edPtEfGIpHVJLeguJrWq2wD4LLApcHAvYjUzsz5o2xTTklqeYjoi5kjag9Sc+nLS0DbnkJJOMc5VCuu+w0tXWwC/zD8PAyYCzwH/Jo1YsBGwCLgZ2C0i7mglTjMz67uWEk6eYvr3pKuEaaTRorcH/hv4sqR9IuLeVo6Zy3c3UgARMaSZdYXti4ADW4nFzMz6j6eYNjOzUrTaLHo4cFK9KaaBk4Gd2hWYmZmtWDzFtJmZlaLVW2rHA2dJejgibu1aWTPF9OfbGZyZvXwNadMoCH0dTWHm6fv3OQ5rD08xbWZmpfAU02ZmVorKp5g2M7OVQ6+GiJW0KTACeBXpVtotEfFYOwMzM7MVS6sdP1cBzgM+zrI9/1+QNAE4KiKWtjE+MzNbQbTaLPoU4KOkxgFDSFNBD8nLH2X5IWnMzMyA1m+p/TdwYmFWz1nANyQFaSDOk9oVnJmZrThavcLZCLirwba78nYzM7PltJpw/g68v8G29wN/61s4Zma2omr1ltpXgEl5sM5LgMdJVzUHA6NonIzMzGwl1+oEbL+Q9Ayp8cC3gFcCi0lTFewTEZPbHqGZma0Qmk44kl5JmnTt7ogYIekVpFk0n3RTaDMz60krdTgvANcB2wJExNKIeMLJxszMmtF0wsmJ5X5g4/4Lx8zMVlSttlL7EnCSpKHtCkDSdpKulbRA0mOSTs0jGnS3z2qSviHpj5IW5j5AjcqOljRd0iJJ90oa067Yzcysea22UjsRWB+4U9KjpFZqy3zYR8Rbmz2YpEHANcC9wGhga+AsUiI8sZtd1wQ+BtwG3ATs3uD4uwK/Ar5N6pS6H/AzSXMi4g/NxmlmZn3XasK5B7i7jec/gjQ8zoERMReYLGkdYLykM/O65UTEM5JeFREh6UgaJBzgy8ANEXF0Xr5e0htJoyE44ZiZlajVZtFj23z+fYGrC4llEnAGsBtweTexNLyNBiBpdVLfoKMLmyYBP5S0bkQ826uozcysZU3V4UgaIOkgScdI+qCkwW06/7bAjNoVETELWJC39cXWpH5CMwrr7yM979f38fhmZtYC9XChgKStSPUsQ2pWzwUO6Ws9iKTFwLERcW5h/Wzgoog4oYljHAmcFxEqrN8FuBHYMSLurFn/WlJru73rxS9pHDAOYPDgwcMmTZrU6tNqm3nz5jFw4MDKzt9If8U1/dG+X3AOHgCPL+zbMYZutu5y6zo1tk6NCzo7tnZY2f4/mzVq1KhpETG83rZmbqmdCSwF3kEaUWBLUiX8d/PvK5SImABMABg+fHiMHDmyslimTJlCledvpL/iGnv87/p8jGOGLuGs6b2aV/BFMw8dudy6To2tU+OCzo6tHVa2/892aOaW2gjSlAR/iohFEXEf8AngNZI26eP55wD1vn4Mytv6emzqHH9QYbuZmZWgmYSzCfBQYd2DgOh7J9AZFOpqJG1OavZcrHtp1YOkcd6KdUHbkq7Y/t7H45uZWQua7fjZfUVP710F7C1p7Zp1Y4CFwNS+HDgingOuJ41kXWsMcLNbqJmZlavZm6NXS1pSZ/21xfUR0cokbBeQmi1fKukMYCvSNNVn1zaVlvQAMDUiDq9Zty+wFrBDXn5f3nR7RDySfz8NmCLpXOAyUsfP/YB9WojRzMzaoJmEc0p/nTwi5kjaAzif1OfmGeAcUtKptSpQHO7mO8AWNcu/zD8PAybm49+YE9FXgE8CDwMf9CgDZmbl6zHhRES/JZx8/HtpPFJAV5khzaxrsO9lpKsbMzOrUKuDd5qZmfWKE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZlcIJx8zMSuGEY2ZmpXDCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccMzMrReUJR9J2kq6VtEDSY5JOlVScTrrefutK+qGkOZKelfQTSesXykyUFHUe2/bfMzIzs3p6nGK6P0kaBFwD3AuMBrYGziIlwhN72P0XwOuBjwFLgTNIU0m/o1BuBnBYYd3MPoRtZma9UGnCAY4ABgAHRsRcYLKkdYDxks7M65YjaQSwF7BbRNyQ1z0K3Cppz4i4pqb4/Ii4pX+fhpmZ9aTqW2r7AlcXEsskUhLarYf9Hu9KNgARcRvwcN5mZmYdpuqEsy3plteLImIWsCBva3q/7L46+20naa6k5yTdKKm7RGZmZv1EEVHdyaXFwLERcW5h/Wzgoog4ocF+k0m3yg4orL8Y2Coi3p6XPwM8T6oj2hA4BhgG7JqviOodexwwDmDw4MHDJk2a1Ovn11fz5s1j4MCBlZ2/kf6Ka/qjz/b5GIMHwOML+3aMoZutu9y6To2tU+OCzo6tHVa2/89mjRo1alpEDK+3reo6nH4VEd+qXZZ0JXAPcAJwQIN9JgATAIYPHx4jR47s3yC7MWXKFKo8fyP9FdfY43/X52McM3QJZ03v29t65qEjl1vXqbF1alzQ2bG1w8r2/9kOVd9SmwPU+/oxKG9r634RsQC4EnhLCzGamVkbVJ1wZlCoc5G0ObAm9etoGu6XNarbqRX5YWZmJao64VwF7C1p7Zp1Y4CFwNQe9ttY0q5dKyQNB7bK2+qSNADYH5jWl6DNzKx1VSecC4DngEsl7Zkr7McDZ9c2lZb0gKTvdy1HxM3AH4CLJB0o6QDgJ8CNXX1w8kgEf5T0CUl7SBoDXA9sCnytpOdnZmZZpY0GImKOpD2A84HLgWeAc0hJp9aqQHG4mzG57A9IifMK4Oia7c8B/yaNWLARsAi4mdRZ9I52Pg8zM+tZ5a3UIuJeYPceygyps+4Z0pA1xWFrurYvAg7se4RmtqIZ0qYWdH1piTfz9P37HMPLTdW31MzMbCXhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZlcIJx8zMSlH5SANWrk7oYQ0rZy9rs5Wdr3DMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSkq7/gpaTvgPGAE8AxwIXBKRLzQw37rAucCB5AS5xXA0RHxVKHcaOArwOuAh/Kxf97WJ2Fm1gYresfsSq9wJA0CrgECGA2cChwDnNLE7r8ARgIfA8YCOwGXFY6/K/Ar4HpgX+B3wM8k7dWO+M3MrHlVX+EcAQwADoyIucBkSesA4yWdmdctR9IIYC9gt4i4Ia97FLhV0p4RcU0u+mXghog4Oi9fL+mNwEnAH/rvaZmZWVHVdTj7AlcXEsskUhLarYf9Hu9KNgARcRvwcN6GpNWBUaQroVqTgBH5lpyZmZWk6oSzLTCjdkVEzAIW5G1N75fdV7Pf1sAr65S7j/S8X9+LeM3MrJcUEdWdXFoMHBsR5xbWzwYuiogTGuw3GZgfEQcU1l8MbBURb5e0C3AjsGNE3FlT5rXA/cDeEbHcbTVJ44BxeXEb4G+9e3ZtsQHwZIXnb6RT4wLH1hudGhc4tt6oOq4tImLDehuqrsPpOBExAZhQdRwAku6IiOFVx1HUqXGBY+uNTo0LHFtvdGpcUP0ttTlAvbqUQXlbX/br+lksN6iw3czMSlB1wplBoa5G0ubAmtSvo2m4X1Zbt/MgsLhOuW2BpcDfexGvmZn1UtUJ5ypgb0lr16wbAywEpvaw38a5nw0AkoYDW+VtRMRzpP43Bxf2HQPcHBHP9j38ftcRt/bq6NS4wLH1RqfGBY6tNzo1rsobDQwC7gXuBs4gJYyzgXMj4sSacg8AUyPi8Jp1V5NGD/g86YrlDOCJiHhHTZldgSnA+aROofvl8vvUazBgZmb9p9IrnIiYA+wBrAJcThph4Bzg5ELRVXOZWmNIV0E/AC4CpgHvLRz/RuB9wJ7A1cB/AR90sjEzK1+lVzhmZrbyqLoOx8zMVhJOOGZmVgonHDMzK4VHGugAkkRq8LA/8AbgVaSWd/8CbgEmRkQl/YZyv6j9AAG/jIinJL2a1Npva2AmMCEippcY0xeAK8s8Z7MkDQBWjYj/1KzbEDgS2I70d70T+PbLpGl+JfL/xHuAt5CmL7mD9DfviErnPKr9k8DuuXFSVTHsDqwG/C4i5uf32qdJLX4fIv1vPlZFfPW40UDF8hvkSmAYKcE8D2xG+ie7ivTG2QY4LSJOKzm2twKTgbWAJcDTwN453heAe4DtgY2BPSPijyXFtZT0+swAfgpMiogHyzh3TyRdCdwfEZ/JyyNIf8elpJaUIv2tnyd9WN1TUlw7AgMi4qaadfsAX+SlRPhXYHxtmZJiuwk4PCLuy8uDSNOHDAPm5WIDSV++9q5N5v0c16e62TwA+AbwLdLYjETEt8uIC14cE/JaYPO86mHSlC2TgfVIHd+3IfVpHBYRs8uKrVsR4UeFD+BnpDfs0Jp1mwK/B36Vl3cj/eN9tOTYJpM6z65HGnn7fGA28BvglbnM6qQP1OtLjGsp8HXSLK/PkZLf7cBngc0q/ns+CYyuWb6F9MGwds26dUlN+q8uMa5bgC/VLH80v47XAl8CTsx/6yW18Zf493xrzfL3SV9u9qlZtw9pOKpzSo7rhfyz3qN22wslv2a/IH1BeC3pjsiP8+fITV3vNdIgnn8FvltmbN3GXXUAK/uDNK32QXXWD8lv6E3y8gnAX0uO7Slg35rljfI/116FcvsDT5YY14sfUPmfbVz+4FySH1PyuldV8PdcALyzZvn54utV85rNLzGuubVxAA8A59Upd0EF77Niwvk38D91yn0eeKTEuH4N/BM4jHw3qGbbejnud5YVT+H8jwGH1CxvkeM5sFDuMODvVcRY7+FGA9V7BSmxFL1Auv3SNfjorVQzh0/U+b14H7ay+7IR8XRETIiIPYBXk6YoX430wflPSX2fJL41d5Mm/uvyOCkpFq1PSk5lWVpY3gK4pE65S0i3Yqq0HqnOpmga6fZtKSLivcBHgGOB2/OUJy9uLiuOBgaRbsF3eTT/fKRQ7iHS/0VHcMKp3mTgK5K26lqR72H/L+kN1dVYYCBQdiXzHcDnJQ2U9ArSVdajwCclrZJjXRX4FOmDtlIR8a+I+FZEvB3YkjRixaYlh3E6cLykj+bX5qvANyS9S9JqklbPdSdfZ/nZaPvTH4FDa5bvAeoNYb8TL314lekgSZ/K9SZzgHrzqWxAulIrTaRRSd5Emin4d5Im5UYzVXuC9KWhywvAd0lfcGptBJRS59WUqi+xVvYH6dvHPaSRrR8gjS23kHSrrfZ21pnAz0uObTjpn39xjukp4M2kJPgQaTiih0n1KKNKjGuZWzCd9gA+RvpgfBa4Lf/+Aul23wv58WtgzRJjGprj+DHwVtJU7E+QEuK7SBXOpwOLqHM7q4S/Z/Hxgzrlvgv8scK/68akYbTmAWflv2NVt9Quq/ca1Sl3HjC5qtes+HArtQ6QrxYOIX2Yr0FKPD+NiKcrDQzI3+beTWpC/6uI+KekjYHjSLdeHgEujIg/lxjTycD3ooOaexZJWp803t9bSR9UIiXv+4ArImJaBTHtAHwH2Jl0S0h5U9fvc4BTI+JbZcfWDEkfBx6MiOsqjmMEaczHbYD9o+RWfTmGwaQvLA/3UO5zpDq5a8uJrHtOOGYrGUlvICWdYiK8KSIWVxmbrdiccDqIpDeSJoirnZV0RpTUV6NVkl4REcXK6MpIWoPUGXUp8EDVH565DmcrajryRsSsKmN6uckdQIkKP6hyZ15FxIKadTuQOz5XcbX6cuVGAx0gVzA/AtwF/JI0gdKE/PtdkmZKOqyi2A6UdJmkKyW9J68bI2kmsFjSI/lWR5kxfUjSR2uWV5V0Oqnvxl2kBgxPSzq+zLhq4hkm6bekytr7gD8BNwMPS3pU0qmS1qwitk4kaa/CJIxIOkDSn0n1h89LukPS/iXHta6kX5PqvuZK+p6kVST9CPgz6f/zNkk3StqgzNiaJekgSfVawVbCCadiko4iVYZeAYwktSp5ZX5sROr0eQVwgaRPlxzbIaRmshuQ/vF/npPLj0n9Xo4mdTS7QNLeJYZ2AqnDaZczcixfB95Jes3OAk6WdEKJcSFpL9JrsinpPv9ppF7zLwDjSRMMHgTclFsjlhnbuyVdK+k+Sb+R9M46ZXau4APqKtKQTl0xvBe4lNSA4fj8WAz8Jr++ZTkNeAfwOVJH2beTWhbuTuqIOphUv7llLms98C21ikl6CLggIs7sodxxwBERsVV35dpJ0u3AtIg4Ii8fSprw7vyIOKam3A+BzSNiz5LiWkBqwTc1Lz8BfLVY2S3p88BREbFFncP0V2zTgLsj4iOF9UeR+ghtReondBNwS0R0N3xKO+N6F2n0iluAvwAjgB2Ac4HPd92ykrQzqS6nOOFhf8a2FHhbRNyWl/8MPBoR7ymUuxJYKyJ2KymumaT31ffy8o6kvkCHRcSPasp9HDghIrYsI658zh80WXQLYGSZf8/u+AqnepuQms725DZK7PSWbcOynQOvIF15FTtTXkoaj6ssz5KuurqsSxrCo+ivpKvEMm0HXFxn/cXAa4BtImIR6YP+vXXK9ZeTgYsiYpeIODIihgEfBz4BXJrrvzrF9qSr/qIJpME8y7IhL/WDgzxmGmmcsloPUL/fUH/6CKkp+9AeHqV92WqGE071/gp8PHesrCtXnH6cVD9RpmDZqb27BlJ8plBuHql3eFl+S+qQulpevgb4QJ1yHyD1cSrTE6Tm7UVvJr2eXZ13H+GlUSTKsD2FRBgRPyDdfnwbcJ2keiMilKX2Vsuz1O+sOJ9yP7MeJr0+Xd5Bavzx9kK5XYCyG4PcD1wXETt19yDdjuwYnp6geseQbnXcK+lS0gjIz+Rt65Jarb2X1EF0n5Jje4T0jf1qgIh4IfdBuK9QbmvSmFNl+SKp5/zdki4kdUA9Q9L2pHHURBpeZkfSEPdlmgCcJmkgKRE+T+q9/yXSAKddfYe2otwPqUWkUb+XERHT8pAtV5Nu840vMaZaV0takn9fl/S3m1oosy3lvs8uAL4laSgpCR5Ceu99Of99/0q64vos5dfh3MLyia+e2v5WlXPCqVhE/Ck3sTyONPTI5oUi/yBVqn4jyh+C/1IK4zBFxK11yr0fKG1OkIh4WtLbSB/in+Ol22Yj8uN50pBB74iI28uKK8f21VwncTxwUtdq0qjg/1NTdDHwtRJDu4s0usBvixsi4qGcdK4EJpYYU5dT6qx7os66g0gjWpciIs7Pdx4+QGoYcFxEXCBpNmnoqa7x8C4AvllWXNl5pJZyPZnKsmP7VcqNBjpMbi67Xl58prbtf6eS9BpSrKWOc1Vz/iEs24nxwQ7og/NK0pXfGsBDVb02NfF8gtS6b8dGI1hIWos05M6eEeHb7d3It7k3iIh/Vx3Ly4kTjpmZlcK31DqE0lTOm5K+nT9ZZ/sGwH4RcVHpwdWR72H/GTi07NtW6vBpnNWB03J3Or1MpkvOVza1U19PI8Vb+jd3ScNJtxlFmoZ+hqQ3k25Rdr3Pzo+Iq8uOrRFf4VRM0uqk1kMH5lVLSSPSfq72w7Ki/hH7dbN5LeDnpLqKuwEi4sqS4urIaZxzLB05LXezlMZZOzgiTi3xnB05XbI6dOrrHMvepMYyT5Na720IjCbVu95H6ms1jNRg5aCIuKys2LpVxpDUfnQ7fPhJpFZpHydNB3A0aU6L+4HX1ZTbmfKnse3IKXbp0Gmc83k7clruFuI/qIL3WUdOl0yHTn2dz/sn0tA6q+TlE3Ic3y+U+zGpg3Hl760ITzFd+YPUDPrIwrqNgRtIU+2OyOuqSDh38NIUu1sUHm/K/5CHdK0rMa6OnMY5n7NTp+V+TZOPIyp4n3XkdMl1Ek5HTH2dz/ks6Qq5a3lQjnf3Qrm9SA16Soutu4frcKq3OYUOnRHxL0l7kL6dXJOHlCmz/0GXnUhXXmeQ+r18PvL8G5K6Oi3+KyKK09r2t65pnG/Iy50yjXOXqPN78d512feyZzZ5TjVZrp1eLtMlr0cHTH2dLWTZflVdvw8olFuT1AerIzjhVO8x4HW89OEJQKRmve+X9C3SpXPpjQUifUWaIOkXwFeA6ZLOy79X6XTgJ5L+QXpduqZxfop0G62r42fZ0zjDS9Ny30hKdrXTcl8XqfNsFdNy/we4Driwh3K7kvqElamTp0s+KFfOQwdNfU26pXaSpPvzub9Jmi34C5JuiIj/5C+Fx5ESYmeo+hJrZX+QBsOc0kOZL1JyPUmDON5E6sn/GPAZqp1it+Omcc5xdeq03JNJQ6H0VK6KOpzL6MDpkungqa9J9V0za97rD5JuiXb9L0wnJec5wA5lxtbdw1c41fs2MEbSq6JBh7yI+LrSfDnvKje05eK4Cxgp6f3AmVQ4ZEZEXJjnKumaxvlpOmAa54i4Iw+FUpyW+528NC33VZQ8LTfpCnpcE+X+TeFquwSfIN366cnDpORUimi+8+sdpBabpYmIB/JQTruQGqdcGxELJY0kfRnbhnRL/qdRUqu+ZrhZtPVKvi20FjAvIjpmgicz61xOOGZmVgqPl/Qykae3/X7VcdTTqbF1alzQ2bF1KknXSLq26jiKOjUu6LzYXIfz8jGKzv2C0KmxdWpc0KGxSbqGdOdjj6pjqUN04GtG58YFHRabb6mZ2Yvyt+FXRETHDGlvK46OyXzWPUlr5GkAOk6nxtapcUHnxhYRe3RqspH0yk58zTo1Lui82JxwXj72JzUL7USdGlunxgUdGltVH1CSPi3pQUkLJf1V0ofrFHsLJb9mnRpXp8fWiBOO2UqiUz+gcr+u80iDsH6Z1IlxoqRLJK1RZiwvh7g6PbbuuNFAxSRd12TRekNq9KtOja1T44LOja3mA+pnpKHr3076gBoNfCgiqhxv6/PANyPixSF18liCPwGul/TuiHjKcb1sYmvIjQYqJmkJ8DfSOEjd2QzYOcqdD6cjY+vUuKBzY5N0B2lom3ofUA8D7440UVwV8y79B3hPREwprB9CGpVhFdI0ABuWGVunxtXpsXXHVzjVuweYERFjuisk6X2UPHwGnRtbp8YFnRvbNqRvxS+KiGslvY30AXWzpH1KjKfWs6QBMJcRETMlvR34HXAzcJrjelEnx9aQ63CqdwvwtibKBeWPXdapsXVqXNC5sTX8gCLdXnuS9AG1U4kxdZkGHFBvQ0TMAfYgjVf2vyXGBJ0bF3R2bA054VTvTOCoJspdCWzZz7EUdWpsnRoXdG5snfwBdTGwlaR6cxoREQuB/yJNrTDLcQGdHVtDrsMxWwlIOhj4LKmupu6o5JJWAb4DvCsiyk7UthJwwjEzs1L4lpqZmZXCCcfMzErhhGNmZqVwwjEzs1I44Vi3JI2VNE3SfyTNkfQXSWf307kOkTS2iXLjJUXN4zFJv5K0dZPnmZh73leu2eecy3Y97/sbbL8/bx/fXzG0eNxlXud2n0fSKyQdmd+TCyXNlXSPpP+V1Ks+TkrulPSRBtsn5t789badL0+q1y0nHGtI0hdJ7fivBg4E/hv4Dal9f384BBjbZNlngRH58XlgB+BaSWs1se9pLZynv7XynAEWAVtKGl67UtJOwJC8vb9jaFbxdW73eX4OfAW4lPSe/Aipf9Pbo/fNbw8BXgX8tBf7fhM4VNJre3nuFZ6HtrHuHAl8NyJOqFl3uaRTqgqoxpKIuCX/foukWcAfgf2AXxYL5z4mq0TE8xHxYIlxttt84M/A+0kdNbu8H7gOGFZFUF3Kep0l7Qu8D9gvIq6q2fTr3l7dZEcDP46IxTXnWpWUPD8MbAp8QNKDwCkR8eLwRHlYmRuBTwLH9CGGFZavcKw76wH/Kq6s/fbYddtE0gGSZkhaJOlGSdsV98u3VKZLek7SPyR9Nf8zI2kicBCwW82tsvEtxDot/xxSJ657SN/8d67dVojtnZKulzRP0rOSpkjasWb7OyRNlbRA0lOSvidp7e4CkjRC0m8l/VPS/Hyr5tDa166Xz3kScEjXB2v+eUhe37YY8mtwSeF4I3OZ7Wtfy55e50bnkbSfpKWStiycZ8u8fnSD12C3/HO50bl7e3WTr0zeDlxS2PQZ4DjSKAxXAh8FfgCsX+cwvyJd5fiztQ5f4Vh3/gwcla8eruhmuPMtgLNJ83IsBE4Brpb0uq5h7yXtRboFchFwLPAm0rfG9YEj8u+vISW5T+Xjzm4h1iH5578K684ETs3r687zImkkMBm4nnRbZj6wC2lE579I2gW4BriM9K16feB0YFBebmQL4E/ABaQP4l2AH0paGhE/o/fP+VLSiAC7kq7q3kEaFfhS4BslxVBrCD2/zo3O80/gMdLrPr6m/FjgCdIglPXMzz+/IemsiHikxZjr2SMf96+F9buRRto+M3+R+lMeg66em4DBwNA6x7GI8MOPug9SUniINNDkUtJIyKcC69SUmZi3v71m3RbAEuCImnW3ANcXjn8c8ALw6rx8CTClibjGkwabXDU/Xk9KFnOBTQpx7VBn/4nAHTXLN5NuT6nB+f5YJ/bd8/G3b/K1VI71u6QPr671TT3n2uedf/8N8H/5928Dl+XfnwTGtyMGYApwSWHdyNrn3eLr3Og8XyElKdXEOZM030uj12Jj4K587gDuBk4ABvbh/T4BuL3O+u8C/8jnnAgM6eYYq+b3/sd7G8eK/PBlnzUUEXcBbyBVyH6b9EHwZeAOSQNrij4RETfV7PcI6RbXW+HF+/pvYfm6lZ+TbuuO6EV46wOL8+NvwFbAmIj4Z02ZRyPizu4OkhsZ7Az8KPInRmH7mjm+X0hatesB3JjP3bDORNIgpRZTj9TEOo6UIPtqEvA+SauTrrKWu51WQgxdenyde/AD0peUkXl5VF7+YaMdIuJfwI7A3qSrvfWArwI3SVoNQNJH8i3EO/Nt3Bn592mSXlnnsBuTEnbRV0lXPg+T/hc+n69668W1BHgmH8sKnHCsWxHxXERcHhFHRsR2wMeA1wGH1xR7os6uTwCb5N83AF4JPF4o07Vcd8TbHjxLGkp/OPBq0rfOqwpliuerZxApkf6zm+2rkBLu4prHc6TntHk3x54IjCHd5torx/sDoB1TAP8WGEj6MFwLuLyCGLo08zo3FBEPka6mDsurDgNui4h7etjvhYj4Q0R8inS77oekW1kj8vYfRcQOpC87S4BdImKHiBgWNY0CaqxB+rsWzzMrH/e9pCv+XYEb1bh7wHO09/VdYbgOx1oSEd+XdCawbc3qjeoU3Yh0Cw7St8bFdcoNzj/rjl7cgyUR0VNfmmYqj+eQbhdu0mD7M/k440kVxkWP1dtJaV75dwOfjogLata35UteRMyXdAVpBOhfRsT8Ypk2xLAIWK2wblC9cJo8XncuBL6n1BT/QFps5RURSyX9gZSsih/2rwPmRM9TLj9NgyuTnKB+rzRV93jSVA/nSDo3J6Ra69G79/QKz1c41pCk5RKJpA2BdVn2W+1GSrMMdpV5Delb5W2QvomSbrEdXDjcIaQP+5vz8vOU/M0wf1DfCvx3V6uvOttvAbaJiDvqPOomHGB10v/Xi9+Yc6u2Yh+mvjzn75CubC5osL2vMcxm2S8WkK6Sequ753pp3j6JFHPdW4QAkgY32PRfwALS37PWm2muAv9v1JmjqN77Arg9/3xVoeyGwJrA35s430rHVzjWnemSfgP8gXSLbAtSJ8sFwI9qyj0JXCzpRF5qpfYE6XZOl5NJLdd+SPowGUpqufS9iOhqFTUDGC3pANKH3WPdfKC30/GkVmhXSZpAul8/glThfQWpccO1kpaSKr7/Q7qFsz/wpYhY7sMlIp6VdDtwkqS5pMR6POlW4Do1RXv9nCPNZz+lm+19jeHXwOGSziG1FhsF9GUa6obPNSIWSfoJ8GngZxHxTDfH+YWk/wC/IDUu2Ag4FBhNqqwv7vtmUgODnvyJ9FptGBH/rln/U0l/AW4g3b4cRrqyfBS4r3CM4aQrvpuw5VXdasGPzn2Q/vn/QLpttIj0z/1TYNuaMhNJLbwOJH2re470j7tc6y1SXcJ00jfZ2aT6h1Vrtm9A+pB7mnwbq0Fc48mttbqJfSI1LaR62kZq+noDKZk+Q2r1tkPN9p2B35Naws0H7iU1BV+3mxheC1yby88iJa5lYm/2ObfwvJdppdbXGIAvklpo/Yc0y+R/sXwrtaZe556eK7BnXr9nD8/xo/lvMTu/l54mJcSRDcpfDry/iff7asBTwIcL69+bz/cvUtKeS0r0O9Y5xrcotGj046WHJ2CzPskd+raPiOE9lTXrTq4bPATYKiKWtvG4s4C9I6J4NVKv7LeA10bE/g22TyQlypl1tq0CPAIcHxEX9ynoFZRvqZlZpSRtA2xHGhLmlDYnm0GkTrHN1ql8A/i7pNdHnVulPTiYdEu5Yf3Tys6NBsysat8l3aq9kjR8TNtExJyIGBCp4Uoz5WeTbtk1arV4GemWaz0CDo/UF8fq8C01MzMrha9wzMysFE44ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NS/D9uOvtcjpd01QAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot probability distribution\n",
    "x = uncertainty_model.values\n",
    "y = uncertainty_model.probabilities\n",
    "plt.bar(x, y, width=0.2)\n",
    "plt.xticks(x, size=15, rotation=90)\n",
    "plt.yticks(size=15)\n",
    "plt.grid()\n",
    "plt.xlabel(\"Spot Price at Maturity $S_T$ (\\$)\", size=15)\n",
    "plt.ylabel(\"Probability ($\\%$)\", size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Payoff Function\n",
    "\n",
    "The payoff function equals zero as long as the spot price at maturity $S_T$ is less than the strike price $K$ and then increases linearly.\n",
    "The implementation uses a comparator, that flips an ancilla qubit from $\\big|0\\rangle$ to $\\big|1\\rangle$ if $S_T \\geq K$, and this ancilla is used to control the linear part of the payoff function.\n",
    "\n",
    "The linear part itself is then approximated as follows.\n",
    "We exploit the fact that $\\sin^2(y + \\pi/4) \\approx y + 1/2$ for small $|y|$.\n",
    "Thus, for a given approximation rescaling factor $c_\\text{approx} \\in [0, 1]$ and $x \\in [0, 1]$ we consider\n",
    "$$ \\sin^2( \\pi/2 * c_\\text{approx} * ( x - 1/2 ) + \\pi/4) \\approx \\pi/2 * c_\\text{approx} * ( x - 1/2 ) + 1/2 $$ for small $c_\\text{approx}$.\n",
    "\n",
    "We can easily construct an operator that acts as \n",
    "$$\\big|x\\rangle \\big|0\\rangle \\mapsto \\big|x\\rangle \\left( \\cos(a*x+b) \\big|0\\rangle + \\sin(a*x+b) \\big|1\\rangle \\right),$$\n",
    "using controlled Y-rotations.\n",
    "\n",
    "Eventually, we are interested in the probability of measuring $\\big|1\\rangle$ in the last qubit, which corresponds to\n",
    "$\\sin^2(a*x+b)$.\n",
    "Together with the approximation above, this allows to approximate the values of interest.\n",
    "The smaller we choose $c_\\text{approx}$, the better the approximation.\n",
    "However, since we are then estimating a property scaled by $c_\\text{approx}$, the number of evaluation qubits $m$ needs to be adjusted accordingly.\n",
    "\n",
    "For more details on the approximation, we refer to:\n",
    "[Quantum Risk Analysis. Woerner, Egger. 2018.](https://www.nature.com/articles/s41534-019-0130-6)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:10.720023Z",
     "start_time": "2020-07-13T23:35:10.711632Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───────┐┌────┐\n",
       "q_0: ┤0      ├┤0   ├\n",
       "     │       ││    │\n",
       "q_1: ┤1 P(X) ├┤1   ├\n",
       "     │       ││    │\n",
       "q_2: ┤2      ├┤2   ├\n",
       "     └───────┘│    │\n",
       "q_3: ─────────┤3 F ├\n",
       "              │    │\n",
       "q_4: ─────────┤4   ├\n",
       "              │    │\n",
       "q_5: ─────────┤5   ├\n",
       "              │    │\n",
       "q_6: ─────────┤6   ├\n",
       "              └────┘</pre>"
      ],
      "text/plain": [
       "     ┌───────┐┌────┐\n",
       "q_0: ┤0      ├┤0   ├\n",
       "     │       ││    │\n",
       "q_1: ┤1 P(X) ├┤1   ├\n",
       "     │       ││    │\n",
       "q_2: ┤2      ├┤2   ├\n",
       "     └───────┘│    │\n",
       "q_3: ─────────┤3 F ├\n",
       "              │    │\n",
       "q_4: ─────────┤4   ├\n",
       "              │    │\n",
       "q_5: ─────────┤5   ├\n",
       "              │    │\n",
       "q_6: ─────────┤6   ├\n",
       "              └────┘"
      ]
     },
     "execution_count": 4,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "# set the strike price (should be within the low and the high value of the uncertainty)\n",
    "strike_price = 1.896\n",
    "\n",
    "# set the approximation scaling for the payoff function\n",
    "c_approx = 0.25\n",
    "\n",
    "# setup piecewise linear objective fcuntion\n",
    "breakpoints = [low, strike_price]\n",
    "slopes = [0, 1]\n",
    "offsets = [0, 0]\n",
    "f_min = 0\n",
    "f_max = high - strike_price\n",
    "european_call_objective = LinearAmplitudeFunction(\n",
    "    num_uncertainty_qubits,\n",
    "    slopes,\n",
    "    offsets,\n",
    "    domain=(low, high),\n",
    "    image=(f_min, f_max),\n",
    "    breakpoints=breakpoints,\n",
    "    rescaling_factor=c_approx,\n",
    ")\n",
    "\n",
    "# construct A operator for QAE for the payoff function by\n",
    "# composing the uncertainty model and the objective\n",
    "num_qubits = european_call_objective.num_qubits\n",
    "european_call = QuantumCircuit(num_qubits)\n",
    "european_call.append(uncertainty_model, range(num_uncertainty_qubits))\n",
    "european_call.append(european_call_objective, range(num_qubits))\n",
    "\n",
    "# draw the circuit\n",
    "european_call.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:11.080005Z",
     "start_time": "2020-07-13T23:35:10.832762Z"
    }
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAE+CAYAAAB1DJw3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA5w0lEQVR4nO3dd5gUVdrG4d8jqCRFFEVcdUbdNeuui2v4TICRYAIVFVcRFXXNcVXWBVQMmHVNKGYUA4qiBAEB0xrANayoqyAgggEFEQaQ8H5/nBrpaXpCz/R0dc2893X1NVPVVd3PNE2/XXXqnCMzwznnnMvGGnEHcM45lzxePJxzzmXNi4dzzrmsefFwzjmXNS8ezjnnsubFwznnXNa8eLiCJqmvJEu5zZY0VNJWMeVpKmmIpB+jPD2i9adJ+krSckkTytm3bdrfUnpbns+/IcqydfTarpe2vkeUqVm+M7lkaRh3AOeq4GfgkOj3LYGrgXGSdjCzRXnOciZwKHAi8A0wVdLGwD3Av4BngHmVPEZ3YFrKchydrbYG+gAPA/NT1r8M7AmU5D+SSxIvHi4JlpvZ29Hvb0uaCbwOdCR8WOfTtsDnZja0dIWkvYEGwINm9lEVHuMjM/tvbQWsCTP7Afgh7hyu8PlpK5dEk6OfxdFppH9J+lxSSXTq6C5J65ZuLOnpTKeSotM230laM1puKemR6JRUiaQJknZN2X46cAqwS8opp76EQgbwYeqprGxFz/ds2rrSU107RsvF0fIxku6T9LOkWZL6SVojbd+dJQ2XNF/SQknvSjpQUltgeLTZV9HjTY/2We20VWWvS+lrI+kmSRdEeeZFp/fWq85r4QqfFw+XRMXRz2+BJoRv/b2BDsCVQHvKHpEMAvaVtEXpCkkCTgIeN7Nl0ephwMHAxUA3wv+P8ZJ+H91/JDAC+IxwamdP4AHgrOj+7tG6lyvJ30BSw5Rbdf4fDgAWAkcBjwP/jH4v/fu2Bd4EWgNnRNmfBzYD3o/+RoAuUeYjK3iuYVT8upQ6Btgf6AX8HegMXFuNv80lgZn5zW8FewP6AnMJp1gbEs7VjwcWAK0zbN8Q2IvQjrB5tG4NYCbQL2W79tE2O0bLh0TL+6Vs05RwCue+lHUPA5PSnrNt6mNV8LeUbpd+uya6fwLwbEWPTSicBjyatt0HwJCU5SeBWUDjcrJ0jh6nOG19j2h9syxfl+nAVKBhyrrbgG/jfg/5rXZu3ubhkmADYFnK8kygm5nNAZD0V+BC4A+ED7ZSWwMzzWylpIeAEyX1tfDJ1oNQBErbHnYDvjeziaU7m9kiSS8Be+f47zmW8EFbanY1HuOVtOUpwOYpy+0JR1WLq/HYqbJ5XcabWeqVY1OAjSStaauO7lwd4cXDJcHPwAGEb8DfArOjAoCkI4FHCVc7XQH8RDhV8zzQKOUxHiKc0mon6T2gK6tO3RDt832G5/4OWD+XfwzwidW8wXx+2vKvlP17NwDm1PA5ILvXJVMmAWtTtvi7OsCLh0uC5WY2qZz7jgbeMbO/la6QtF/6RmY2XdJYwhHHFoRTWU+mbDIH2CjD47ciFKR8WAKslbauRTUf60fCB39NFcLr4gqQN5i7pGsMLE1b172cbQcRjjj+Bgwzs/kp971DOMWyb+kKSU2ATsAbOUtbsVmES4FTHVTNxxoHHCOpUTn3/xr9LO/+UoXwurgC5MXDJd0YwpVUvSUdIOkWwhU/mQwjfLv/M+E01m/MbDTwFvCUpJMkdSZcWdUYuLG2wqd5HviDpFujv6U/qzpHZqsf0Bx4TVK36PEukdQzuv/z6OfpknaXtFOmBymQ18UVIC8eLunuA24GzgOeA4qA4zNtaGZLgZHA18DYDJscQShGtxEu9RXQ3sy+zHXocvK9TGi3OYpQSIoIf1d1HutzQoP2XMLlxM9Hjzsjun8Goc2nC+GS3uGZHwmI+XVxhUlRu6NzdZ6khoQPzwfN7Mq48ziXZN5g7uo8SWsBfyQckWxAOFpxztWAFw9XH2wCvEu45PR0M5sVcx7nEs9PWznnnMuaN5g755zLWr04bdWyZUsrLi6u1r6LFi2iadOmlW9YIJKUN0lZIVl5k5QVkpU3SVmhZnknT54818w2zHhn3INr5ePWpk0bq67x48dXe984JClvkrKaJStvkrKaJStvkrKa1SwvaYOApt78tJVzzrmsefFwzjmXNS8ezjnnsubFwznnXNa8eDjnnMuaFw/nnKuLBg+G4mL2a98eiovDcg7Vi34ezjlXrwweDL16QUkJApgxIywDdC9vupvs+JGHc87VNb17Q0lJ2XUlJWF9jnjxcM65umbmzOzWV4MXD+ecq2uaN8+8fvPNc/YUXjycc64ueeopmD8fGjQou75JE+jfP2dP48XDOefqitdegxNPhH32gQcegKIiTIKiIhg4MGeN5eBXWznnXN0wZQocfjhsuSUMGwbrrw89ejBxwgTatm2b86fzIw/nnEu6OXOgQwdo1AhGjgyFo5b5kYdzziXZL79Ap07w44/htFU15y7KlhcP55xLqmXL4Oij4aOPYPhw+POf8/bUXjyccy6JzOCMM2D06NA43qFDXp/e2zyccy6Jrr4aHnwQrrwSTjkl70/vxcM555LmoYegTx846STo1y+WCF48nHMuSV55JQxyeOCBcP/9IMUSw4uHc84lxQcfQNeusMMO8OyzsOaasUXx4uGcc0kwcyZ07AgtWsDLL8O668Yax6+2cs65QjdvXriaqqQE3nwTfve7uBN58XDOuYK2dCkceSR88UW4LHeHHeJOBHjxcM65wrVyJfToARMnhtkB27WLO9FvvM3DOecK1eWXw5AhcP31cPzxcacpw4uHc84VorvuggED4Mwz4dJL406zGi8ezjlXaF54Ac49Fw47DO68M7a+HBXJe/GQtL2kcZJKJM2WdJWkBlXYb1dJr0j6KbqNlbR7PjI751zevPMOHHcc7LorPPnk6jMCFoi8Fg9JLYCxgAGHA1cBFwEV9q+XtFm0X0Pgr9GtITBGUlFtZnbOubz58kvo3Blatw6j5DZpEneicuX7aqszgMZAFzNbQPjwXxfoK2lAtC6TTsA6wJFm9jOApLeAuUBH4J7aj+6cc7Xohx9CXw4zGDUKNtoo7kQVyvdpqw7A6LQiMYRQUParYL81geXAopR1C6N1hXcy0DnnslFSAoceCrNmwUsvwR/+EHeiSuW7eGwLfJa6wsxmAiXRfeUZGm1zs6SNJG0E3ArMA56ppazOOVf7VqwIl+G++y488QTssUfciaok38WjBTA/w/p50X0ZmdlsoB3QFfguunUBDjazH3If0znn8sAMzjsvXF11xx2hJ3lCyMzy92TSMuASM7stbf0s4FEzu6Kc/VoDrwFTWNW+cRawC/B/0dFL+j69gF4ArVq1ajNkyJBqZV64cCHNmjWr1r5xSFLeJGWFZOVNUlZIVt5cZt1syBC2uu8+ZnbrxrQzzsjJY6arSd527dpNNrNdM95pZnm7Ad8DfTKsX0QoKuXtdwswHVgzZd1awAzgjsqet02bNlZd48ePr/a+cUhS3iRlNUtW3iRlNUtW3pxlffJJMzDr1s1sxYrcPGYGNckLTLJyPlfzfdrqM9LaNqLLcJuQ1haSZlvgEzNbVrrCzH4FPgG2qoWczjlXeyZODLMA7rsvPPwwrJG8/tr5TjwSOFjSOinrugGLgYkV7DcD2FHSWqUrJK0N7Eg4InHOuWSYMgWOOAK22gqGDYNGjeJOVC35Lh73AkuB5yQdELVL9AVusZTLdyV9KWlQyn4PAJsAz0vqJKkzMAxoDQzMV3jnnKuR2bNDX45GjWDkyDCxU0LltXiY2Txgf6ABMJzQs/xWoE/apg2jbUr3mwwcQugo+BjwKOFU14Fm9mHtJ3fOuRr65Rfo1Al++glGjICiZA+Okff5PMxsCtC+km2KM6wbB4yrpVjOOVd7li2Do46Cjz8OU8juskvciWrMJ4NyzrnaZAannw6vvAKDBsHBB8edKCeS18TvnHNJctVV8NBD0KcP9OwZd5qc8eLhnHO15cEHoW9fOPnkUDzqEC8ezjlXG0aPhl694KCD4L77CnJCp5rw4uGcc7n2n/+EBvKddoJnnoE114w7Uc558XDOuVyaMQM6doT11w9XVq27btyJaoVfbeWcc7kyb17oBLh4MYwdC5tsEneiWuPFwznncmHp0jDsyNSpob1jhx3iTlSrvHg451xNrVwZBjp87TV48klo2zbuRLXO2zycc66mLrsMnnoKbrgBjj027jR54cXDOedq4q674MYb4ayz4JJL4k6TN148nHOuuoYNg3POgcMPh9tvr3N9OSrixcM556rj7bfhuONgt93giSegQYPK96lDvHg451y2vvwSDj0UNt0Uhg+HJk3iTpR3Xjycc64qBg+G4mL2a98ettsOliwJEzptuGHcyWLhxcM55yozeHAYp2rGDGQGy5eH2zvvxJ0sNl48nHOuMr17Q0lJ2XVLloT19ZQXD+ecq8zMmdmtrwe8eDjnXGWaN8+8fvPN85ujgHjxcM65ijz5JMyfv/qluE2aQP/+sUQqBF48nHOuPBMmQI8esN9+Yf7xoiJMgqIiGDgQunePO2FsfGBE55zL5JNPwii5v/89PP88tGgBJ53ExAkTaFsPBj6sjB95OOdcutmzw7wcTZrAiBGhcLgy/MjDOedSLVgQZgKcNw9efz2conKr8eLhnHOlli0Lc4//979hCtk//SnuRAXLi4dzzgGYhV7kY8bAgw/CwQfHnaigeZuHc84B9OsHDz8MffvCySfHnabgefFwzrlBg0Lx6NkT/vnPuNMkghcP51z9NmoUnH56OE117731akKnmvDi4Zyrv95/H44+GnbeGZ55BtZcM+5EiZH34iFpe0njJJVImi3pKklVmoJLUhdJ70laLOlHSaMkNa3tzM65Omj6dOjUCdZfP1xZtc46cSdKlLxebSWpBTAWmAIcDmwF3EwoYv+oZN9TgX8BA4BLgBZAe/yKMedctubNC305liyBceOgdeu4EyVOvj94zwAaA13MbAEwRtK6QF9JA6J1q5HUErgVOMfM7k+56/laT+ycq1uWLAnDjkydGi7L3X77uBMlUr5PW3UARqcViSGEgrJfBfsdE/18pLaCOefqgZUr4aST4LXX4JFHYN99406UWPkuHtsCn6WuMLOZQEl0X3l2Bz4HTpE0S9IySe9I+r/ai+qcq3P+/nd4+mm48UY49ti40yRavotHC2B+hvXzovvKszGwDaFd5O/AocAiYJSkVjnO6Jyri+68E266Cc4+Gy66KO40iSczy9+TScuAS8zstrT1s4BHzeyKcvZ7BTgQ6GBmo6J16wIzgH+Z2ZUZ9ukF9AJo1apVmyFDhlQr88KFC2nWrFm19o1DkvImKSskK2+SskLt5235+uvs0KcPc/fai0/69l19Yqcs1KfXtl27dpPNbNeMd5pZ3m7A90CfDOsXEYpKefs9BawEGqWtHwsMrex527RpY9U1fvz4au8bhyTlTVJWs2TlTVJWs1rO+9ZbZo0ame2xh9miRTV+uPr02gKTrJzP1XyftvqMtLYNSZsBTUhrC0nzKaDoVmZ3QlFxzrnV/e9/cOihsOmm8OKLYX4OlxP5Lh4jgYMlpfbG6QYsBiZWsN9L0c92pSskNQfaAB/mOqRzrg74/vswoZMEI0fChhvGnahOyXfxuBdYCjwn6YCoXaIvcIulXL4r6UtJg0qXzWwS8AIwSNJJkjoBLwLLgLvy+Qc45xJg0aJwxDFnDrz0UphK1uVUXouHmc0D9gcaAMOBfoTOf33SNm0YbZPqBGAYcAvwLKFwtI8e0znnguXL4bjjYNIkGDIEdt897kR1Ut6H9jCzKYRhRSrapjjDuoXAmdHNOedWZwbnngvDh8Ndd8Fhh8WdqM7yUXWdc3XHgAFwzz2hM+Df/hZ3mjrNi4dzrm544gm47LJwyuraa+NOU+d58XDOJd/48dCjB7RtCw89BGv4R1ttq/QVlvSqpG2j30+UtEHtx3LOuSr673/hyCNh663h+edh7bXjTlQvVKU87wOsF/3+EGEODueci98334R5OZo2hREjYL314k5Ub1TlaquvgaMlLST06N4i+j2j6Goq55yrXQsWhMIxbx68/jpsvnncieqVqhSP64C7gfMBA54oZztF91d/xDHnnKuKZcvgqKNgypQwheyf/hR3onqn0uJhZvdLehH4A/AacBZhGlnnnMs/MzjttDAL4EMPwUEHxZ2oXqq0eEg6EXjZzN6Q1A94wcxm134055zLoE+fMAtgv37hCisXi6o0mKc2kv8T2LT24jjnXAUeeACuvhpOOQWuXG0aH5dHVSke84BNot9L2zWccy6/Ro6EM86AQw4JvciVPkODy6eqNJiPBR6T9Hm0/LCkReVtbGa75SSZc86VmjwZjj4adt45zEG+5ppxJ6r3qlI8ehIGI9wW+DPwFfBDbYZyzrnfTJ8OnTpBy5bhyqp11ql0F1f7qnK1VQlwM4CkA4DeZuYTMDnnat9PP4UJnX79NQxB0rp13IlcJKsh2c1si9oK4pxzZSxZAocfDtOmwdixsN12cSdyKbIePUzSlpLukfSxpG+in3dL2rI2Ajrn6qGVK+Gkk+CNN+Cxx2CffeJO5NJkdeQhqQ0wHlhCmFf8O6AV0BXoLqmdmb2f85TOufrl0ktDw/hNN8Exx8SdxmWQ7UyCNwH/ATpEbSEASGoCjIjur3CWQOecq9Add8DNN8M558CFF8adxpUj29NWuwEDUgsH/NaofhPgkwU756rvuefg/PPDEOu33up9OQpYtsVjMVDefB7rE05nOedc9t56C7p3hz32gMGDoYGPsVrIsi0eLwPXS9o7dWW0fB0wPFfBnHP1wODBUFzMfu3ahUbx5s3hxRehceO4k7lKZFs8LgSmARMlzZH0oaQ5wERC58GLch3QOVdHDR4MvXrBjBkIwhVWP/8Mo0fHncxVQbb9PH4E9pZ0CPAXoDUwB3jHzF6phXzOubqqd28oKSm7bsmSsL5793gyuSrL9lLdBma2wsxGAaNqKZNzrj6YOTO79a6gZHva6htJAyR5V0/nXPWZhXnHM/HpZBMh2+JxL3AU8F9J70jqJWndWsjlnKvLbrgBFi6EhmknP5o0gf7948nkspJV8TCzvma2JXAg8DlwCzBH0uBo0ETnnKvY4MFw+eVw/PFhGtmiIkyCoiIYONDbOxIi67GtAMzsVTM7EdgYOAfYBhgtabqkvpI2qfgRnHP10quvwsknQ7t28OCDcMIJMH06E199NQy97oUjMapVPFLsCuxLmOtjHvA6cCrwpaQTavjYzrm65OOPQ8/xrbcOPcnXXjvuRK4GqjOqbpGkPpKmAuMIl+v2BDYxs78CRcB9wI05TeqcS65Zs6BjR2jWLEwnu956cSdyNZTtpbrjgX2Ab4CHgIfMbEbqNma2QtITwHk5S+mcS64FC8JMgD//DK+/DpttFncilwPZHnl8D3QEiqPG8xnlbPcBkHHiKEnbSxonqUTSbElXSaryIDaS1pA0SZJJ6pxlfudcPv36K3TtClOmwNCh8Mc/xp3I5Ui2Pcy7VXG7ZcBqhUVSC2AsMAU4HNiKMMXtGsA/qhjjVGDTKm7rnIuLGZx2WpgF8OGH4cAD407kcijb+TwAkLQpsDXQKP0+MxtRwa5nAI2BLma2ABgT9RPpK2lAtK6i520B9AcuAx6oTnbnXJ7885/w6KNw9dVhVkBXp2Tb5rEO8DRwUOmq6KelbFbRKagOwOi0IjEEuAHYj8pH5b0aeJPQUO+cK1T33w/XXAOnnhrGqnJ1TrZtHtcBmxMazQUcCbQFBhFG1d2jkv23BT5LXWFmM4GS6L5ySdqZcFXXxVlmds7l04gRcOaZ0KED3HOPT+hUR8nMKt+qdGNpGqFt4ilgGbC7mb0X3XczsJmZlTvhsKRlwCVmdlva+lnAo2Z2RQX7TiSM3nuppGJCsTrUzF4qZ/teQC+AVq1atRkyZEiV/85UCxcupFmzZtXaNw5JypukrJCsvHFlbfb55+xy/vmUbLYZH9x+OyuqOC+Hv7a1pyZ527VrN9nMds14p5lV+QYsAvaJfv8FODjlvv2B+ZXsvww4P8P6WcC1Fex3LPAtsG60XEw4Vda5KrnbtGlj1TV+/Phq7xuHJOVNUlazZOWNJeu0aWatWpkVF5vNmZPVrv7a1p6a5AUmWTmfq9metvoaaBn9/gWQeqns7lQ+De08oHmG9S2i+1YjaU1Ch8MbgDUkrQeUDsbYNGqHcc7F6ccfw2mqX38NnQA33jjuRK6WZXu11RjgAOB54FbgEUltgKWEYUpurmT/z0hr25C0GdCEtLaQFE0Jl+beEt1SDQGmAr+v+p/gnMupJUvg8MPD2FRjx8K2FTZfujqiSsVDUmNC58DvgWmSWpnZY5IWEoZobwycTRiWpCIjgUskrWNmv0TrugGLCVPZZrIQaJe2bmPgSeAK4NWq/A3OuVqwciX89a/w1lvw1FOw995xJ3J5UmnxkLQloWNfccrqBZKOMbPnCUchVXUvcC7wnKQbgC2BvsAtlnL5rqQvgYlmdoqZLQcmpGUqzfKxmb2TxfM753Lp4ovh2Wfh5pvh6KPjTuPyqCptHgOAlYTLc5sAOwD/ofKjjNWY2TxCw3oDQp+OfoTTX33SNm1Ixf1FnHNxu/12uPVWOO88uOCCuNO4PKvKaas9gYvM7M1o+VNJp0c/W5vZnGye0MymAO0r2aa4kvuns6qDonMu34YODQWjS5dw1OF9Oeqdqhx5tAampa2bSvjw9ksqnKtv3nwzTNq0557w+OPQwE8S1EdVvVS36j0JnXN11+efw2GHhSljX3gBqtgJ0NU9Vb1Ud7Sk5RnWj0tfb2Yb1TyWc67gfPdd6MvRsGHoy9GyZeX7uDqrKsWjX62ncM4VtkWLoHPnUEAmTIAtt4w7kYtZpcXDzLx4OFefLV8O3brB+++HU1V/+UvciVwBqNZ8Hs65esIMzjoLXn4Z7r03HH04R/ZDsjvn6pPrroOBA+Hyy+H00+NO4wqIFw/nXGaPPx4mcjrhBOjfP+40rsB48XDOrW7cOOjZE9q3h0GDvBOgW40XD+dcWR9/HHqOb7MNPPccrLVW3IlcAfLi4ZxbZdas0JdjnXXCdLLNM02/45xfbeWcK/Xzz9CxI/zyC7z+Omy2WdyJXAHz4uGcCzMAdu0Kn34Ko0bBzjvHncgVOC8eztV3ZnDqqaGR/JFHYP/9407kEsDbPJyr7668Eh57DK65Bk48Me40LiG8eDhXnw0cGPpwnHYaXHFF3GlcgnjxcK6+evllOPPM0Eh+993el8NlxYuHc/XRpElwzDGwyy7w1FNhmHXnsuDFw7n6Zto06NQJNtoIXnoJmjWLO5FLIC8eztUnP/4YOgEuWxYmdNrYZ5J21ePHqs7VF4sXhylkZ8yAsWNh223jTuQSzIuHc/XBihXw17/Cv/8NTz8Ne+8ddyKXcF48nKsPLr4Yhg6FW26Bo46KO42rA7zNw7m67rbbwu388+GCC2IO4+oKLx7O1WXPPgsXXhjGrbr55rjTuDrET1s5V5cMHgy9e7PfzJnhUty5c2HPPcPwI2v4d0WXO148nKsrBg+GXr2gpAQBfPdd6DV+4onQuHHc6Vwd419FnKsreveGkpKy68zguuviyePqNC8eztUVM2dmt965GvDi4VxdUd7Mf5tvnt8crl7Ie/GQtL2kcZJKJM2WdJWkBpXs8xdJD0n6Mtrvc0l9JDXKV27nCpoZbLXV6uubNAlDrjuXY3ktHpJaAGMBAw4HrgIuAvpVsms3YCvgBqAjcBdwITC41sI6lyTXXgvjx4fhR4qKMAmKisJ8Hd27x53O1UH5vtrqDKAx0MXMFgBjJK0L9JU0IFqXyfVmNjdleYKkJcB9korMbEYt53aucD32GPzjH2H4kUceAYmJEybQtm3buJO5Oizfp606AKPTisQQQkHZr7yd0gpHqf9EPzfJXTznEmbsWOjZM8w7/sADPqGTy5t8F49tgc9SV5jZTKAkui8bewIrgam5ieZcwnz0EXTpAtttF8atWmutuBO5ekRmlr8nk5YBl5jZbWnrZwGPmlmVJlGWtDHwETDCzHqUs00voBdAq1at2gwZMqRamRcuXEizBE2Wk6S8ScoKhZV37e+/589nnQXA+3ffzdINNyxzfyFlrYok5U1SVqhZ3nbt2k02s10z3mlmebsBy4DzM6yfBVxbxcdYC3gNmAa0qMo+bdq0seoaP358tfeNQ5LyJimrWQHlnTfPbMcdzdZd1+yjjzJuUjBZqyhJeZOU1axmeYFJVs7nar4bzOcBzTOsbxHdVyFJAh4FdgD2MrNK93GuTvn113Cq6vPPw0yAO+0UdyJXT+W7eHxGWtuGpM2AJqS1hZTjNsIlvgeaWVW2d67uMAuN4+PHw6OPhkZy52KS7wbzkcDBktZJWdcNWAxMrGhHSZcDZwMnmNkbtRfRuQLVu3cY/LB//3BZrnMxynfxuBdYCjwn6YCoUbsvcIulXL4b9SQflLJ8PHAt4ZTVN5L2SLmVbSl0ri66774wwGGvXnD55XGncS6/p63MbJ6k/YF/AcOB+cCthAKSnit1yJKDop89oluqk4GHcxrUuULy0kvwt79Bp05w113el8MVhLzP52FmU4D2lWxTnLbcg9WLhnN133vvQbdu8Oc/w1NPQUOfgscVBh9V17lCNW0adO4MrVqFo4+mTeNO5NxvvHg4V4jmzoVDDoHly8Mlua1axZ3IuTL8GNi5QrN4cRgdd+ZMGDcOttkm7kTOrcaLh3OFZMUKOOEEePtteOYZ2GuvuBM5l5EXD+cKhRlceCE89xzcdht07Rp3IufK5W0ezhWKW2+FO+6ACy6A886LO41zFfLi4VwheOYZuOgiOOoouOmmuNM4VykvHs7F7fXXw3Aje+8dZgVcw/9busLn71Ln4vTpp3D44VBcDC+8AI0axZ3IuSrx4uFcXL79Fjp0CDMAjhwJ668fdyLnqsyvtnIuDgsXhrGq5s6FCRNgiy3iTuRcVrx4OJdvy5fDMcfAhx/Ciy/Crpln+XSukHnxcC6fzODMM8NpqoEDoWPHuBM5Vy3e5uFcPvXvDw88AP/4B5x2WtxpnKs2Lx7O5csjj8CVV8KJJ8JVV8Wdxrka8eLhXD6MGQOnngoHHAD33+8TOrnE8+LhXG378MMwTtX228Ozz4ZLc51LOC8eztWmr78OjeLNm8OIEeGnc3WAX23lXG2ZPz90Aly4EN54A373u7gTOZczXjycqw1Ll8KRR8L//gejRsFOO8WdyLmc8uLhXK6tXAk9e4ae448/Du3bx53IuZzzNg/ncq13b3jiCbj2WujePe40ztUKLx7O5dI998D118MZZ8Bll8Wdxrla48XDuVx58UU4+2zo3BnuvNP7crg6zYuHc7nw7rtw7LHQpg0MGQINvTnR1W1ePJyrqalTw9FG69YwfDg0bRp3IudqnRcP52pi7tzQl2PlyjBSbqtWcSdyLi/82Nq56lq8GA47LPQiHzcOtt467kTO5Y0XD+eqY8WKcBnu22+H8ar+7//iTuRcXnnxcC5bZnDBBfD883D77dClS9yJnMu7vLd5SNpe0jhJJZJmS7pKUoMq7Ndc0kOS5kn6WdJgSRvkI7NzZdxyS7gU98IL4dxz407jXCzyeuQhqQUwFpgCHA5sBdxMKGL/qGT3p4GtgVOBlcANwDBgn1qK69zqnnoKLr44zEF+441xp3EuNvk+8jgDaAx0MbMxZnYv0A+4UNK65e0kaU/gIOAkMxtqZs8DJwB7SzqgVpIOHgzFxezXvj0UF4flQpakvEnKCqvytmsX+nJss02YFXANv1jR1V/5fvd3AEab2YKUdUMIBWW/Svb7zsxeK11hZu8CX0X35dbgwdCrF8yYgcxgxoywXKgfcknKm6SsUDZv6bqZM2Ho0DhTORe7fDeYbwu8mrrCzGZKKonuG17Bfp9lWP9pdF9u9e4NJSVl15WUQI8eYbC7QvO//8Hy5WXXFWreirJefXVYNiv7M9O6qv6s6T4//BD6cKRavDi8R3zQQ1eP5bt4tADmZ1g/L7qvOvttmWkHSb2AXgCtWrViwoQJVQ6538yZZBqVyJYv54cNN6zy4+TLhlOmJCZvhVlbt/5tPCjLNC5U6brytkm7vyrbVPY8rYcPz5x35kwmZvGeyreFCxdm9Z6PW5LyJikr1GJeM8vbDVgGnJ9h/Szg2gr2GwMMy7D+ceCtyp63TZs2lpWiIrPw3bPsragou8fJlyTlTVJWs+TljYwfPz7uCFlJUt4kZTWrWV5gkpXzuZrvNo95QKZJnFtE9+V6v+rp3x+aNCm7rkmTsL4QJSlvkrJC8vI6lyf5Lh6fkdZGIWkzoAmZ2zTK3S9SXltIzXTvDgMHQlFROK1RVBSWC/Ucd5LyJikrJC+vc3mS7+IxEjhY0jop67oBi4GJley3saS9S1dI2pXQ3jGyNoLSvTtMn87EV1+F6dML/8MiSXmTlBWSl9e5PMh38bgXWAo8J+mAqFG7L3CLpVy+K+lLSYNKl83s38ArwKOSukg6AhgMvGFmY/P5BzjnnMtz8TCzecD+QAPCZbn9gFuBPmmbNoy2SdWNcHTyIPAoMBk4sjbzOuecyyzvAyOa2RSgfSXbFGdYNx84Obo555yLkY+v4JxzLmtePJxzzmVNljokQx0l6QdgRjV3bwnMzWGc2pakvEnKCsnKm6SskKy8ScoKNctbZGYZh6moF8WjJiRNMrNd485RVUnKm6SskKy8ScoKycqbpKxQe3n9tJVzzrmsefFwzjmXNS8elRsYd4AsJSlvkrJCsvImKSskK2+SskIt5fU2D+ecc1nzIw/nnHNZ8+LhnHMua148nHPOZc2Lh3POuax58XDOOZe1vI+q63IjmoGxIyDgGTP7UdKmwMXAVsB0YKCZfRxfSpD0d2BE3DmqSlJjoKGZ/ZKybkPgbGB7YCXwAXC3mf0cS0jnCoBfqhuRJML8IJ2A7YD1CR8U3wJvAw+b2f/iS7iKpN2AMUBTYDnwE3AwMAJYAXwC7AhsDBxgZq/HFBVJKwEjTBf8BDDEzKbGlacykkYAX5jZedHynoTZKlcS5pAR0Ab4FWhvZp/EmHUXoLGZvZWy7hDgclYVug+BvqnbFIro/9yhwJ8J75FJhC8aBf2hJGldwlhR7c3sjbjzwG+Z2gNrAS+b2aLoS89ZhBlXpxG+TM7O2XMW+L9TXkQv8gjCh8K3hA+G3xHe0CMJL/42wNVmdnVcOUtJGkM4ajwSWESYUOsIwofbUWa2TNLawDCgkZm1iylqafG4AdgJOJCQ+31CIXnazL6JK1smkuYCp5jZC9Hy24TX+IjSoxFJzYEXgSVmdnCMWd8GhptZ/2i5J/AAMB54lVDo9gf2AbqW/k0xZX2L8Lp+Gi23IMwO2gZYGG3WjPBF7eDUI784SPpbBXc3Bm4Ebge+ADCzu/ORKxNJvwfGAZtFq74CDiJ8wVwPmEr4/FoMtDGzWTl5YjOr9zfgScKbYKeUdZsAo4Ch0fJ+hDd5zwLI+yPQIWV5I8K3zIPStusEzI0560pgt+j39YFe0Rt9eXSbEK1bP+7XNcpYAuybsvxr+uua8touijnrgtRswJfAnRm2uxf4sFDeB9HyIMIR8yEp6w4B5gG3FsD7YCXhKH5lObfU+1bEnPVpwhHm76P/Y49Fn2dvAetE27SMtrkvV8/rDeZBB+AySzkvb+Hw7gzgCEmtzWwicC1wXkwZ01mG39MPIwvqsNLMfjKzgWa2P7ApcBHhMPteYI6kl2MNGPwXSD1S+47wHzLdBoRCE6eVactFwLMZtnuW8M2zkBwGXGVmo0pXRL/3B7rElmqVF4DvgVOABma2RumN8H4Q0DZalz5ldr7tDfQ3sy/N7CfgH4R2z5ssOoIzs7nAbZR9b9eIF49gDcI3iXQrCG+S5tHyO8DW+QpVgUnAxZKaSVoDuAL4BjhTUgMASQ2BvxE+DAuOmX1rZreb2f8BWxDmsd8k5lgA1wOXSeoZvYb9gRslHShpLUlrR+0K1xG+8cXpdaB7yvInQKaht/9CeH8UkvUI7+N0kwltdbEysyOBk4BLgPck7ZV6dzypytWCcLq9VOm/dfocRtMIX9pywq+2CsYA10j6yMymwW/nZO8g/KOUNpQ3AwrhCpvehMzzCKd+SgiNZc8AX0gqbTDfhHAqoKCZ2QzCh/b1BZDlOUnnEL6l3Qp8TvjyMIrwoaFo0xcJHyxxugJ4M/oCcSehofwRSesTTgeK8L44H7gspoypukoqLW7zgEyTDLUknI6LnZm9Imlnwuv3sqRRhKsZY22PyeB7wlFnqRXAfYSj5lQbkcPs3mAORJe4jiYcVcwgnOfeAlgKHGdmI6PtBhBm1uoWV9ZSUebOhC8AQ81sjqSNgUsJpyhmAA+Y2fsxxkRSH+B+y+FVHvkgaQOgG7Ab4ZuwCB94nwIvmdnkGOP9RtKfgHuA3Slb3Ep/n0c4PXR7LAEj0YUT6R42s55p290HbG9m++QnWdVE/7cGEE6p3UcoKO3M7LU4cwFIGgb8lP5aZtjuTmBbMzswJ8/rxSOITvccA/wRaERofHwiOofoXEGTtB2hgKQXurfMbFmc2bIh6TRgqpm9GneWTKJLt28lfEHrZAVwCbSkVkATM/uqku0uJFw4MS4nz+vFo+6RtIaZZfqmVzAkNSI06q0EvizED7iozWNLUvr8mNnMeFM5Vxi8wTyNpB0kdZV0anTrKmmHuHOlk9RF0jBJIyQdGq3rJmk6sEzSjOhbXKwknRD1PyhdbijpesJlmh8RGvR/klQI5+QBkNRG0ouE88OfAm8C/wa+kvSNpKskNYk1ZB2iSNw5MpHUOP3fWtKfos+FNnHlKghxXp9cSDegJ6GdINO13SsIw32cHHfOKOsxUa43CJcUlgCnEdpqBhF6lT4Z5T445qxTgDNTlm+O8l4J7EW4zLAvoQPTFQXw2h5EaOuaRLg0uy+ho+ivUeaLCFc1fQC0KIC8nQn9Zj6N3gv7Zthmd+Lvi3AQUZ+DlHVHEDqMLgeWRa95p7hf0yhbc+D5KNdy4H6gAfBI2ufCG0DLuPNW8W/qmsv3Qex/UCHcgHOiN8ldhN64LaM3SoPo972Bf0UfKmcVQN73gHtTlrtH2W5O2+4hYGzMWUuA/VKWvwfOy7DdxcCMAnhtJwOPlPMemU44Wm8UfejdHXPWA6MPsDej9+fkaPlmolPS0XaFUDxWULaT4JHRB/Bb0b/9RYSju+Vk6JQZQ947CEOQnAOcGH1hGAp8HRXCDQn9w74B7ok7bxX/ppwWD2/zACRNI3wYD6hku0uBM8xsy/wkKzfHAqCLmY2NlpsTGkgPsJSGxuh01n1mFlv/CUlzgLPNbGi0vJRwNDQhbbsDgRfNrHH+U5bJsRg4zMzGpK1vQejZv4OZfSrpROAGM2sdR84o0xuEcbhOTlnXk/DBN4ZwpeASSbsTGs5j68wWXW21h5m9Gy2/D3xjZoembTcCaGpm+8UQMzXHdELHu/uj5V0IxflkM3skZbvTCEfMW8QSNGR4sIqbFhE6NubkfeBtHkFr4N0qbPcuBdCBiXAZZuoboHRsoPlp2y0kdMaK04uEDo1rRctjgeMybHcc4dtd3L4nXHGX7o+E1720n88MVnUejcuOwOOpK8zsQcJQOnsAr0Z9PgrRjoRLXtMNJAyUGLcNWdW/C6IxrAjjRKX6ksz9VfLpJMLR0E6V3IrKe4Dq8E6CwYfAaZJes3KuUooa9E4jNPLGbQZh1NTRAGa2IrqE8NO07bYC5uQ5W7rLCT2h/yvpAWA4cIOkHVnVka0dsAthhNW4DQSultSMUOh+JfTQ7g2Mt1X9VbYE4r7yaglhZOUyzGxy1CN6NOG0UN885ypP6mmOn8ncYW0RhfGl9itCEZ4YLe9DOM32f4R2jlJ7Ef/74AvgXTM7saKNJB0FPJWrJ/XiEVxE6EE8RdJzhOHD50f3NQe2JZyj3ZTC6LH9HGnDDJjZOxm2O5ayb/S8M7OfJO1B+PC9kNDLFWDP6PYr4RTLPmb2XjwpVzGz/tEplsuAf5auJlyAcH7KpssIDepx+ohw3v3F9DvMbFpUQEYAD+c5V3lGS1oe/d6c8IVhYto22xL/Fx4I463dLmknQqE7hvBF6Mroi8WHhCOkC4C4R9p+m1DUKpPaibTGvM0jImkrQu/sQ1g1tHGprwlX3NxoBTwXRTpJmwPzzawghnsAkFRM2Y5sU60w+3isSThyawRMK6TXsJSk0wlDlOxi5XRmldSUcNXQARYG9YtFNNJAui/M7Im07SZE6wvhMvNzCadT1ySM1nCvpOMIbUqlA2MOBP4e53s4umR4LzO7o5LtWhLa7NILdvWe14vH6qLruteLFuebWdyjpzrnCkR0Crulmf0Qd5Y4efGoY6JD6veB7oVwGkgJnNZVCZni17k4efFIEX1obEI4lTI3w/0tgY5m9mjew5XN0bGCu5sSGsUuIxqO3cxG5CNXJkrQtK6QrCl+qyoa9+poM7sq5hx5nyo1l6IjjtRpcycT/o7YP0QVRivuSvj/9LCZfSbpj0A/Vn3h+ZeZjc7Zk8bdcaUQbsDahOHMV0S3ZYSe2s3Ttou9s1WUI0mznM0FDk9ZfpvQI3qdlHXNCQ2nowvgtR1DmMZ1PcK57n8Bswi9t9dMeb+MJFx9Ffv7twp/U047h1Uzw+8JVwmWvi+nEj7UphEK9HuEodi/AzYtgNfsLWC7lOUWUcaVUc4FrOrkuE5cOaNsBxO+fH0bva4LCFcwzovy3RX9v1tBmE45N88b9z9SIdwIV9XMJ1yKuytwbvQm/gL4Q8p2hVI8JhGuSDmZcO126m3n6E19TOm6mLMmZlrXKEeSpvjdvIq3M+J+3xLTVKk1yJuYaXMJIww8Q5jxEMJFFPOAQWnbPQa8nbPnjfsfqRBuhEtzz05btzHwGvADsGe0rlCKhwjzfn9PGDJhi5T7mkdv/NXGOIop67tAn5Tlr4FjM2x3IvBDAeT9Me0DYsPo9TwwbbuOBVA8So8yK7sVwhHobOCYlOWiKFeXtO1OBv5XAO+D9OLxA3B+hu1iH1aHcCnxASnLLaL87dO2O4hwAVBOntf7eQSbkdb5z8y+lbQ/oVqPldSdwrj+HAvvhIGSngauAT6OJnq5Jt5kGV0PDJb0NfAoq6Z1/ZFwqqq0k2AhTOsKq6b4fYNw1JQ6xe+rFjpkFsoUv78ArwIPVLLd3oTL0OMUy1SpObQehTtt7mLKdhYt/T19qJ8mhI6lOeHFI5gN/IFwpPEbC9duHyvpdsJhYawN5enMbD5wtqSBhGvPvwBuoIDmWLZkTesKyZri911Cu9zLFW0UzZ0St1imSq2hpEyb+ybwT0lfRFluIoxm/fdo1IxfovHvLiUUu5zwq634bWCxLc2sbQXbXE741mwW4wBzFZF0LGGqzE0JA6DFPkVmKSVkWldI1BS/VwK9zCy9U2v6dvsC/cysXX6SZcwwjBimSq0uJWjaXEm/JwylU/o+mE44mn+W0GN/BlBM+DLUzsw+yMnzevH47TK3bsB1VsG0s5KOJ5z7Prm8beIWnVJpCiw0sxVx53EO4psqtbapQKbNjfp37UW4QnCcmS2OOjufyqovPE+Y2aycPacXD+ecc9kqhNErXS2RdL+kQXHnqIokZYXk5XUu17zBPAuS7gfWMLNT4s5SRe1IzheEJGWFBOWVNJZwlmH/uLNUJklZIVl5c53Vi0d2EvOBAWBmv487Q1UlKSskLq9Izvs2SVkhWXlzmtXbPOqw6BLNjcws7slqKpWkrJC8vM7lWlIqZkGQ1CiaIyMpOhFmREuCJGWFBOWVtGZS3rdJygrJypvrrF48spOYDwxXP0g6S9JUSYslfSjprxk2+zMF8L5NUlZIVt44snqbRwJJquo15Zl6xOZVkrJCsvJGnULvJEyR+x/CVKQPSzocOMHMcjYURU0lKSskK29cWb14kKwPjMi+hGE+plSyXSEMS5GkrJCsvBcDN5nZb+NWReOxDQbGS+psZj/Glq6sJGWFZOWNJas3mAOSllO1D4zfAbvHPTyJpA+Bz8ysWyXbHQU8FWfeJGWNciQmr6RfgEPNbELa+mLCfCMNCONvbQi85VmrLkl548rqbR7BJ8B/zezoim7ALXEHjbwN7FGF7VIHHoxLkrJCsvL+TBiYrwwzm044dTEX+Dfwl/zGyihJWSFZeWPJ6kce/Da42SFmVlTJdl0Jc1rHWnQlbQXsYGYvVrJdY8LlpOnDXudNkrJGORKTV9ILwC9mdkI59zcmDI7XgZgH9ExS1ihPYvLGldWLB8n6wHCulKSjgQuAzuUN6CmpAXAPYUDPLfKZLy1HYrJGWRKTN66sXjycc85lzds8nHPOZc2Lh3POuax58XD1jqQekiZL+kXSPEn/kVQrV9JJ2lpSX0nrVWHbvpIs5TZb0tCoTa6yfR+WlGmObedqhRcPV68oTCf8ADAa6AKcCLwAHFZLT7k10AdYr4rb/wzsGd0uBv4EjJPUtJL9rgZ6VCuhc9XgPcxdfXM2cJ+ZXZGybrikfnEFSrPczN6Ofn9b0kzgdaAj8Ez6xpIam9liM5uaz5DO+ZGHq2/WA75NX2kplx1KKo5OGx0v6bHo9Nb3kvqk7yepvaR3JC2R9J2kuxXmk0ZSW2B4tOlX0WNOzzLv5OhncfSY0yXdLOlKSbOABdH61U5bSSqS9KSkuZJKJH0k6fiU+xtJGiDpa0lLowH1OmaZz9VTfuTh6pv3gXOib/QvVTLmz43AS8BRhDGv+kiaa2Z3AUjaARgFjAG6ApsB1wNbEoaDeJ9o3CHCKbI5wNIs8xZHP1ML3vGEURH+Rjn/hyVtROhVXBJl+BrYMcpY6llgN8JptanAMcCLknY1sw+yzOnqGS8err45CxgGPAyYpE+BoYSB5RakbfuJmZ0e/T46+kC+QtI9ZrYSuBKYARxmZisAJP0EPCVpTzP7t6TPo/3/Ew0XUSlJpf8vtwTuBn4BxqZt1rmS0VIvAJoDbcxsTrRuXMpz7E+YYqCtmU2MVr8iaWugN3B0VbK6+stPW7l6xcw+ArYjNJDfTRif6kpgUunpphTPpy0/B2wCbBot7wY8X1o4IkOB5cDe1Yy4AbAsun1OKCDdUgoAwLgqDLPdHhiVtl+qAwhHM29Kalh6IxSYXauZ3dUjfuTh6h0zW0poixgOIOkUwhVYpwC3p2z6fdqupcutgZnRz+/SHnuFpB+B9asZ72fCB7sRPtxn2+rDQHy32l6r2wB4r4L7WwIbE4pUuhUZ1jlXhhcPV++Z2SBJA4Bt0+7aqJzlOSk/y2wTjSG0AZBxjKEqWG5mlfXXqMqYQj8Silt5fgK+AY6oYi7nyvDTVq5eidot0tdtSGgfSP9Gf2Tacmmj96xo+R3gyKhgpG7TEHgjWv41+pnvyaPGAQdLalXB/RsDC81sUvotfzFdUvmRh6tvPo6GsH6FcBqqiHA1UgnwSNq2O0TD9Q8lXG11CnBe1FgOcA1h2s9hku4htIXcAIw2s39H25Q2mJ8uaQhQYmYf186fVsathA6Qr0vqT7jaajugqZkNIFwhNhoYI+kGwtVb6xI6JTYys8vzkNElmBcPV99cBRwO3EFol/gWeIvQKP1V2raXAp0JxWMJoRf3v0rvNLNPJHUAriU0pi8gzCN9aco2MyRdDJwLnEM4aimujT8slZn9IGkvYABwG7A28AVwXXS/SeoCXAGcD2xOOJX1AWE+bOcq5EOyO5cmmr7zK8LUni/FHMe5guRtHs4557LmxcM551zW/LSVc865rPmRh3POuax58XDOOZc1Lx7OOeey5sXDOedc1rx4OOecy9r/AymBnmomK8ydAAAAAElFTkSuQmCC\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "# plot exact payoff function (evaluated on the grid of the uncertainty model)\n",
    "x = uncertainty_model.values\n",
    "y = np.maximum(0, x - strike_price)\n",
    "plt.plot(x, y, \"ro-\")\n",
    "plt.grid()\n",
    "plt.title(\"Payoff Function\", size=15)\n",
    "plt.xlabel(\"Spot Price\", size=15)\n",
    "plt.ylabel(\"Payoff\", size=15)\n",
    "plt.xticks(x, size=15, rotation=90)\n",
    "plt.yticks(size=15)\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:11.087435Z",
     "start_time": "2020-07-13T23:35:11.082555Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact expected value:\t0.1623\n",
      "exact delta value:   \t0.8098\n"
     ]
    }
   ],
   "source": [
    "# evaluate exact expected value (normalized to the [0, 1] interval)\n",
    "exact_value = np.dot(uncertainty_model.probabilities, y)\n",
    "exact_delta = sum(uncertainty_model.probabilities[x >= strike_price])\n",
    "print(\"exact expected value:\\t%.4f\" % exact_value)\n",
    "print(\"exact delta value:   \\t%.4f\" % exact_delta)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate Expected Payoff"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───────┐┌────┐\n",
       "q_0: ┤0      ├┤0   ├\n",
       "     │       ││    │\n",
       "q_1: ┤1 P(X) ├┤1   ├\n",
       "     │       ││    │\n",
       "q_2: ┤2      ├┤2   ├\n",
       "     └───────┘│    │\n",
       "q_3: ─────────┤3 F ├\n",
       "              │    │\n",
       "q_4: ─────────┤4   ├\n",
       "              │    │\n",
       "q_5: ─────────┤5   ├\n",
       "              │    │\n",
       "q_6: ─────────┤6   ├\n",
       "              └────┘</pre>"
      ],
      "text/plain": [
       "     ┌───────┐┌────┐\n",
       "q_0: ┤0      ├┤0   ├\n",
       "     │       ││    │\n",
       "q_1: ┤1 P(X) ├┤1   ├\n",
       "     │       ││    │\n",
       "q_2: ┤2      ├┤2   ├\n",
       "     └───────┘│    │\n",
       "q_3: ─────────┤3 F ├\n",
       "              │    │\n",
       "q_4: ─────────┤4   ├\n",
       "              │    │\n",
       "q_5: ─────────┤5   ├\n",
       "              │    │\n",
       "q_6: ─────────┤6   ├\n",
       "              └────┘"
      ]
     },
     "execution_count": 7,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "european_call.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:11.253873Z",
     "start_time": "2020-07-13T23:35:11.247765Z"
    }
   },
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "problem = EstimationProblem(\n",
    "    state_preparation=european_call,\n",
    "    objective_qubits=[3],\n",
    "    post_processing=european_call_objective.post_processing,\n",
    ")\n",
    "# construct amplitude estimation\n",
    "ae = IterativeAmplitudeEstimation(\n",
    "    epsilon_target=epsilon, alpha=alpha, sampler=Sampler(run_options={\"shots\": 100, \"seed\": 75})\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:35.384882Z",
     "start_time": "2020-07-13T23:35:11.373858Z"
    }
   },
   "outputs": [],
   "source": [
    "result = ae.estimate(problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:35.392375Z",
     "start_time": "2020-07-13T23:35:35.387058Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:        \t0.1623\n",
      "Estimated value:    \t0.1687\n",
      "Confidence interval:\t[0.1637, 0.1737]\n"
     ]
    }
   ],
   "source": [
    "conf_int = np.array(result.confidence_interval_processed)\n",
    "print(\"Exact value:        \\t%.4f\" % exact_value)\n",
    "print(\"Estimated value:    \\t%.4f\" % (result.estimation_processed))\n",
    "print(\"Confidence interval:\\t[%.4f, %.4f]\" % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Instead of constructing these circuits manually, the Qiskit Finance module offers the `EuropeanCallPricing` circuit, which already implements this functionality as a building block."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit_finance.applications.estimation import EuropeanCallPricing\n",
    "\n",
    "european_call_pricing = EuropeanCallPricing(\n",
    "    num_state_qubits=num_uncertainty_qubits,\n",
    "    strike_price=strike_price,\n",
    "    rescaling_factor=c_approx,\n",
    "    bounds=(low, high),\n",
    "    uncertainty_model=uncertainty_model,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:        \t0.1623\n",
      "Estimated value:    \t0.1687\n",
      "Confidence interval:\t[0.1637, 0.1737]\n"
     ]
    }
   ],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "problem = european_call_pricing.to_estimation_problem()\n",
    "# construct amplitude estimation\n",
    "ae = IterativeAmplitudeEstimation(\n",
    "    epsilon_target=epsilon, alpha=alpha, sampler=Sampler(run_options={\"shots\": 100, \"seed\": 75})\n",
    ")\n",
    "result = ae.estimate(problem)\n",
    "\n",
    "conf_int = np.array(result.confidence_interval_processed)\n",
    "print(\"Exact value:        \\t%.4f\" % exact_value)\n",
    "print(\"Estimated value:    \\t%.4f\" % (european_call_pricing.interpret(result)))\n",
    "print(\"Confidence interval:\\t[%.4f, %.4f]\" % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Evaluate Delta\n",
    "\n",
    "The Delta is a bit simpler to evaluate than the expected payoff.\n",
    "Similarly to the expected payoff, we use a comparator circuit and an ancilla qubit to identify the cases where $S_T > K$.\n",
    "However, since we are only interested in the probability of this condition being true, we can directly use this ancilla qubit as the objective qubit in amplitude estimation without any further approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:35.826269Z",
     "start_time": "2020-07-13T23:35:35.819320Z"
    }
   },
   "outputs": [],
   "source": [
    "from qiskit_finance.applications.estimation import EuropeanCallDelta\n",
    "\n",
    "european_call_delta = EuropeanCallDelta(\n",
    "    num_state_qubits=num_uncertainty_qubits,\n",
    "    strike_price=strike_price,\n",
    "    bounds=(low, high),\n",
    "    uncertainty_model=uncertainty_model,\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">         ┌──────┐\n",
       "state_0: ┤0     ├\n",
       "         │      │\n",
       "state_1: ┤1     ├\n",
       "         │      │\n",
       "state_2: ┤2     ├\n",
       "         │  cmp │\n",
       "state_3: ┤3     ├\n",
       "         │      │\n",
       " work_0: ┤4     ├\n",
       "         │      │\n",
       " work_1: ┤5     ├\n",
       "         └──────┘</pre>"
      ],
      "text/plain": [
       "         ┌──────┐\n",
       "state_0: ┤0     ├\n",
       "         │      │\n",
       "state_1: ┤1     ├\n",
       "         │      │\n",
       "state_2: ┤2     ├\n",
       "         │  cmp │\n",
       "state_3: ┤3     ├\n",
       "         │      │\n",
       " work_0: ┤4     ├\n",
       "         │      │\n",
       " work_1: ┤5     ├\n",
       "         └──────┘"
      ]
     },
     "execution_count": 14,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "european_call_delta._objective.decompose().draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<pre style=\"word-wrap: normal;white-space: pre;background: #fff0;line-height: 1.1;font-family: &quot;Courier New&quot;,Courier,monospace\">     ┌───────┐┌──────┐\n",
       "q_0: ┤0      ├┤0     ├\n",
       "     │       ││      │\n",
       "q_1: ┤1 P(X) ├┤1     ├\n",
       "     │       ││      │\n",
       "q_2: ┤2      ├┤2     ├\n",
       "     └───────┘│  ECD │\n",
       "q_3: ─────────┤3     ├\n",
       "              │      │\n",
       "q_4: ─────────┤4     ├\n",
       "              │      │\n",
       "q_5: ─────────┤5     ├\n",
       "              └──────┘</pre>"
      ],
      "text/plain": [
       "     ┌───────┐┌──────┐\n",
       "q_0: ┤0      ├┤0     ├\n",
       "     │       ││      │\n",
       "q_1: ┤1 P(X) ├┤1     ├\n",
       "     │       ││      │\n",
       "q_2: ┤2      ├┤2     ├\n",
       "     └───────┘│  ECD │\n",
       "q_3: ─────────┤3     ├\n",
       "              │      │\n",
       "q_4: ─────────┤4     ├\n",
       "              │      │\n",
       "q_5: ─────────┤5     ├\n",
       "              └──────┘"
      ]
     },
     "execution_count": 15,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "european_call_delta_circ = QuantumCircuit(european_call_delta._objective.num_qubits)\n",
    "european_call_delta_circ.append(uncertainty_model, range(num_uncertainty_qubits))\n",
    "european_call_delta_circ.append(\n",
    "    european_call_delta._objective, range(european_call_delta._objective.num_qubits)\n",
    ")\n",
    "\n",
    "european_call_delta_circ.draw()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:35.843101Z",
     "start_time": "2020-07-13T23:35:35.828358Z"
    }
   },
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "problem = european_call_delta.to_estimation_problem()\n",
    "\n",
    "# construct amplitude estimation\n",
    "ae_delta = IterativeAmplitudeEstimation(\n",
    "    epsilon_target=epsilon, alpha=alpha, sampler=Sampler(run_options={\"shots\": 100, \"seed\": 75})\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:58.524643Z",
     "start_time": "2020-07-13T23:35:35.848588Z"
    }
   },
   "outputs": [],
   "source": [
    "result_delta = ae_delta.estimate(problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:58.531941Z",
     "start_time": "2020-07-13T23:35:58.526967Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact delta:    \t0.8098\n",
      "Estimated value: \t0.8091\n",
      "Confidence interval: \t[0.8034, 0.8148]\n"
     ]
    }
   ],
   "source": [
    "conf_int = np.array(result_delta.confidence_interval_processed)\n",
    "print(\"Exact delta:    \\t%.4f\" % exact_delta)\n",
    "print(\"Estimated value: \\t%.4f\" % european_call_delta.interpret(result_delta))\n",
    "print(\"Confidence interval: \\t[%.4f, %.4f]\" % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:35:58.967675Z",
     "start_time": "2020-07-13T23:35:58.792732Z"
    }
   },
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Software</th><th>Version</th></tr><tr><td><code>qiskit</code></td><td>None</td></tr><tr><td><code>qiskit-terra</code></td><td>0.45.0.dev0+c626be7</td></tr><tr><td><code>qiskit_finance</code></td><td>0.4.0</td></tr><tr><td><code>qiskit_aer</code></td><td>0.12.0</td></tr><tr><td><code>qiskit_algorithms</code></td><td>0.2.0</td></tr><tr><td><code>qiskit_ibm_provider</code></td><td>0.6.1</td></tr><tr><td><code>qiskit_optimization</code></td><td>0.6.0</td></tr><tr><th colspan='2'>System information</th></tr><tr><td>Python version</td><td>3.9.7</td></tr><tr><td>Python compiler</td><td>GCC 7.5.0</td></tr><tr><td>Python build</td><td>default, Sep 16 2021 13:09:58</td></tr><tr><td>OS</td><td>Linux</td></tr><tr><td>CPUs</td><td>2</td></tr><tr><td>Memory (Gb)</td><td>5.778430938720703</td></tr><tr><td colspan='2'>Fri Aug 18 16:00:58 2023 EDT</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    },
    {
     "data": {
      "text/html": [
       "<div style='width: 100%; background-color:#d5d9e0;padding-left: 10px; padding-bottom: 10px; padding-right: 10px; padding-top: 5px'><h3>This code is a part of Qiskit</h3><p>&copy; Copyright IBM 2017, 2023.</p><p>This code is licensed under the Apache License, Version 2.0. You may<br>obtain a copy of this license in the LICENSE.txt file in the root directory<br> of this source tree or at http://www.apache.org/licenses/LICENSE-2.0.<p>Any modifications or derivative works of this code must retain this<br>copyright notice, and modified files need to carry a notice indicating<br>that they have been altered from the originals.</p></div>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import tutorial_magics\n",
    "\n",
    "%qiskit_version_table\n",
    "%qiskit_copyright"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "celltoolbar": "Tags",
  "kernelspec": {
   "display_name": "Python 3",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.9.7"
  },
  "varInspector": {
   "cols": {
    "lenName": 16,
    "lenType": 16,
    "lenVar": 40
   },
   "kernels_config": {
    "python": {
     "delete_cmd_postfix": "",
     "delete_cmd_prefix": "del ",
     "library": "var_list.py",
     "varRefreshCmd": "print(var_dic_list())"
    },
    "r": {
     "delete_cmd_postfix": ") ",
     "delete_cmd_prefix": "rm(",
     "library": "var_list.r",
     "varRefreshCmd": "cat(var_dic_list()) "
    }
   },
   "types_to_exclude": [
    "module",
    "function",
    "builtin_function_or_method",
    "instance",
    "_Feature"
   ],
   "window_display": false
  },
  "vscode": {
   "interpreter": {
    "hash": "e3b168dd14084693aa742087410f9921d6040e41eb6bdb17b20e4003862f82dd"
   }
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}