{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pricing Bull Spreads"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "<br>\n",
    "Suppose a [bull spread](http://www.theoptionsguide.com/bull-call-spread.aspx) with strike prices $K_1 < K_2$ 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",
    "\n",
    "$$\\min\\{\\max\\{S_T - K_1, 0\\}, K_2 - K_1\\}$$\n",
    "\n",
    "\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",
    "\n",
    "$$\\mathbb{E}\\left[ \\min\\{\\max\\{S_T - K_1, 0\\}, K_2 - K_1\\} \\right]$$\n",
    "\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",
    "$$\n",
    "\\Delta = \\mathbb{P}\\left[K_1 \\leq S \\leq K_2\\right]\n",
    "$$\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://arxiv.org/abs/1806.06893)\n",
    "- [Option Pricing using Quantum Computers. Stamatopoulos et al. 2019.](https://arxiv.org/abs/1905.02666)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "import matplotlib.pyplot as plt\n",
    "\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\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": {},
   "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 circuit for uncertainty model\n",
    "uncertainty_model = LogNormalDistribution(\n",
    "    num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high)\n",
    ")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "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_1$, then increases linearly, and is bounded by $K_2$.\n",
    "The implementation uses two comparators, that flip an ancilla qubit each from $\\big|0\\rangle$ to $\\big|1\\rangle$ if $S_T \\geq K_1$ and $S_T \\leq K_2$, and these ancillas are 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",
    "\n",
    "$$ \\sin^2( \\pi/2 * c_\\text{approx} * ( x - 1/2 ) + \\pi/4) \\approx \\pi/2 * c_\\text{approx} * ( x - 1/2 ) + 1/2 $$\n",
    "\n",
    "for small $c_\\text{approx}$.\n",
    "\n",
    "We can easily construct an operator that acts as\n",
    "\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",
    "\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://arxiv.org/abs/1806.06893)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set the strike price (should be within the low and the high value of the uncertainty)\n",
    "strike_price_1 = 1.438\n",
    "strike_price_2 = 2.584\n",
    "\n",
    "# set the approximation scaling for the payoff function\n",
    "rescaling_factor = 0.25\n",
    "\n",
    "# setup piecewise linear objective fcuntion\n",
    "breakpoints = [low, strike_price_1, strike_price_2]\n",
    "slopes = [0, 1, 0]\n",
    "offsets = [0, 0, strike_price_2 - strike_price_1]\n",
    "f_min = 0\n",
    "f_max = strike_price_2 - strike_price_1\n",
    "bull_spread_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=rescaling_factor,\n",
    ")\n",
    "\n",
    "# construct A operator for QAE for the payoff function by\n",
    "# composing the uncertainty model and the objective\n",
    "bull_spread = bull_spread_objective.compose(uncertainty_model, front=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAE+CAYAAAB1DJw3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA7tklEQVR4nO3dd5xU1f3/8ddb0QiiiA0riyX+bKmQGEuk2MEGqKgYJRqxt68lUTRgwdh7RbElGOy9IB27AVuioqICgiiiIOIiUj6/P84ZGcbZ3Zlldu7c3c/z8ZjH7r1z78x7h2E+c++55xyZGc4551wxVkg6gHPOufTx4uGcc65oXjycc84VzYuHc865onnxcM45VzQvHs4554rmxcNVNEkDJFnW7TNJD0naLKE8q0oaKumrmKdPXH+0pE8kLZI0poZ9O+X8LZnbonL+DTHLFvG1XSNnfZ+YqWW5M7l0aZZ0AOcK8A2wZ/x9U+BCYKSkbczsuzJnOQ7YBzgcmA58JGk94GbgBuABYHYdj9Eb+DhrOYnOVlsA/YG7gDlZ658Ctgeqyx/JpYkXD5cGi8zslfj7K5KmAs8DXQkf1uW0JfC+mT2UWSFpJ2BF4A4ze7uAx3jbzP7XUAGXh5l9CXyZdA5X+fy0lUujCfFnu3ga6QZJ70uqjqeObpS0emZjSffnO5UUT9t8IWmluLy2pLvjKalqSWMkdcjafjJwFPCbrFNOAwiFDOCt7FNZxYrP92DOusyprm3jcru4fJCkWyV9I2mapPMlrZCz7y8lPSFpjqR5kl6TtJukTsATcbNP4uNNjvv85LRVXa9L5rWRdIWk02Ke2fH03hr1eS1c5fPi4dKoXfz5OdCC8K2/H7AXcB7QhWWPSAYDO0vaJLNCkoAjgH+Z2cK4+lFgD+AMoBfh/8doSZvH+7sDTwMTCad2tgduB06I9/eO656qI/+Kkppl3erz//AyYB5wAPAv4O/x98zftyXwIrA+cGzM/giwMfB6/BsBesTM3Wt5rkep/XXJOAjYBegL/BXYG7i4Hn+bSwMz85vfKvYGDABmEU6xNiOcqx8NzAXWz7N9M2BHQjtC27huBWAqcH7Wdl3iNtvG5T3jcsesbVYlnMK5NWvdXcD4nOfslP1Ytfwtme1ybxfF+8cAD9b22ITCacA9Odu9CQzNWv43MA1oXkOWvePjtMtZ3yeub1nk6zIZ+AholrXuGuDzpN9DfmuYm7d5uDRYC1iYtTwV6GVmMwAk/Qn4P+DnhA+2jC2AqWa2RNKdwOGSBlj4ZOtDKAKZtoffAzPNbGxmZzP7TtKTwE4l/nsOJnzQZnxWj8d4Lmf5XaBt1nIXwlHV/Ho8drZiXpfRZpZ95di7wLqSVrKlR3eukfDi4dLgG2BXwjfgz4HPYgFAUnfgHsLVTucAXxNO1TwCrJL1GHcSTml1lvQfoCdLT90Q95mZ57m/ANYs5R8DvGPL32A+J2f5B5b9e9cCZiznc0Bxr0u+TAJ+xrLF3zUCXjxcGiwys/E13Hcg8KqZHZ9ZIalj7kZmNlnSCMIRxyaEU1n/ztpkBrBunsdvQyhI5fA9sHLOutb1fKyvCB/8y6sSXhdXgbzB3KVdc2BBzrreNWw7mHDEcTzwqJnNybrvVcIplp0zKyS1ALoBL5Qsbe2mES4FzrZ7PR9rJHCQpFVquP+H+LOm+zMq4XVxFciLh0u74YQrqfpJ2lXSVYQrfvJ5lPDt/reE01g/MrNhwEvAfZKOkLQ34cqq5sDlDRU+xyPAzyVdHf+WgSztHFms84FWwDhJveLjnSnpyHj/+/HnMZK2k/SLfA9SIa+Lq0BePFza3QpcCZwCPAxUAYfm29DMFgDPAJ8CI/Jssj+hGF1DuNRXQBczm1Tq0DXke4rQbnMAoZBUEf6u+jzW+4QG7VmEy4kfiY87Jd4/hdDm04NwSe8T+R8JSPh1cZVJsd3RuUZPUjPCh+cdZnZe0nmcS7OyH3lI2jz2jH1b0uJ8PX/z7PM7SXdKmhR7uL4vqX8t53Od+5GklSX9jnCaZS3C0YpzbjkkcbXVNoQxiV4BVipwn17AZsClwIfALwmD4/2S0ADqXG02AF4jXHJ6jJlNSziPc6lX9tNWklYwsyXx9weBtc2sUx37rG1ms3LW9SV8g2wXz98655wrk7KftsoUjiL3mZVn9Rvx5wbLl8g551yx0txJcHtgCcsO85DX2muvbe3atavXk3z33XesuuqqdW9YIdKUN01ZIV1505QV0pU3TVlh+fJOmDBhlpmtk+++VBYPhcl3zgX+aWb5hk7InNbqC9CmTRuuuOKKej3XvHnzaNkyPZOqpSlvmrJCuvKmKSukK2+assLy5e3cuXPNTQJJjsoIPAiMKXKflYFxhJnYWheyT/v27a2+Ro8eXe99k5CmvGnKapauvGnKapauvGnKarZ8eckZQTr7lqojjzgHwz2EK7Z2NLO6pvt0zjnXAFJVPAg9XPcDdjOziQlncc65Jis1xUPS2cCJwEFm5gOyOedcgspePOKInF3j4obA6pIy02c+bWbVkiYBY83sqLjPoYTpLO8Cpkv6Q9ZDfmRmX5YnvXPOOUhmYMR1CYOrPQD8Adg6azkzb0AzwrzUGZlhqfsAL+fcujV4YuecGzIE2rWjY5cu0K5dWK5kDZy37EceZjaZMCpnbdu0y1nuQygczjlXfkOGQN++UF0dPrymTAnLAL1rmj4mQWXIm5o2D+ecS0y/flBdvey66mo44QR4//38+yTpuuvy5+3Xz4uHc86VzdSp+dd/8w1cdFF5sxSipjELa/o76sEng3LOubqstVb+9VVVsGRJ5d2qqvLnbdu2ZC+JFw/nnKvNtdfCrFmwQs7HZYsWMHBgMpnqMnBgyJetxHm9eDjnXD5mcM45cOqp0KMHDB4MVVWYFL7ZDxpUmY3lEHINGtSgeb3Nwznnci1aBMccA3fcEX7eeCOsuCL06cPYMWPo1KlT0gnr1rs39O7dYHn9yMM557LNnw89e4bC8fe/w803h8LhluFHHs45lzF7Nuy7L7z4YjjaOP74pBNVLC8ezjkHMH067LknfPAB3HcfHHhg0okqmhcP55x7/33Yffdw5PHMM9ClS9KJKp4XD+dc0/baa9C1a2jXGDMGfvvbpBOlgjeYO+earueeC0cZrVqFdg4vHAXz4uGca5ruvRe6dYPNNw+FY/PNk06UKl48nHNNz7XXhn4QO+0EY8fCeuslnSh1vHg455oOMzj77KW9xp95JpyyckXzBnPnXNNQU69xVy9+5OGca/yqq8ORxh13QP/+3mu8BPzIwznXuM2eDfvsAy+95L3GS8iLh3Ou8fJe4w3Gi4dzrnHyXuMNyouHc67x8V7jDa7sDeaSNpd0q6S3JS2WNKbA/VpJulPSbEnfSBoiqYa5IZ1zTdawYd5rvAySuNpqG6Ar8D7wQRH73Q90Av4C9AF+Bzxa2mjOuVS7917Ye2/vNV4GSZy2esLMHgOQ9CCwdl07SNoe2B3oaGbj4rrpwKuSdjWzEQ0Z2DmXAtdeGzr/deoEjz7qnf8aWNmPPMxsST122wv4IlM44uO8BnwS73PONVXeazwRaWkw3xKYmGf9e/E+51xT5L3GEyMzS+7J42krM+tUx3bDge/MbP+c9f8CNjWzHfLs0xfoC9CmTZv2Q4cOrVfGefPm0bJly3rtm4Q05U1TVkhX3jRlhfrlXeH779n6wgtZ+6WXmHzEEUw+4giQGijhUk3htc3o3LnzBDPrkPdOM0vsBjwIjClgu+HAo3nW/wt4qa7927dvb/U1evToeu+bhDTlTVNWs3TlTVNWs3rk/fprsx13NJPMbryxQTLVpNG/tlmA8VbD52paTlvNBtbJs751vM8511RMnw577AEffgj33w8HHJB0oiYpLQMjTiR/20ZNbSHOucZo4kTYYQeYOjU0jHvhSExaisczwHqSdsqskNQB2DTe55xr7F57LUze9P33ode4DzeSqLKftpLUgtBJEGBDYHVJma8PT5tZtaRJwFgzOwrAzF6W9Bxwj6QzgCXApcAL5n08nGv8hg2Dnj2hTZvwu3f+S1wSbR7rAg/krMssbwJMJuTKvd6uF3A1cAfhiOlJ4OQGS+mcqwxDhkCfPrDttuFUlU8ZWxHKXjzMbDJQ6/V0ZtYuz7o5wJ/jzTnXFFxzDZx2mvcar0BpafNwzjUlmV7jp50WTld5r/GKk5ZLdZ1zTcWiRdC3L9x5Jxx7LNxwg/car0B+5OGcqxyZucbvvDPMNX7TTV44KpQfeTjnKkP2XOM33QTHHZd0IlcLLx7OueQMGQL9+tFx6lRo1gyWLPFe4ynhxcM5l4whQ0LbRnV1uPxy4UL42c9gwYKkk7kCeJuHcy4Z/fqFNo5sCxaE9a7iefFwziVj6tTi1ruK4sXDOVd+Q4aEvhz5tG1b3iyuXrx4OOfK65pr4LDDYKutoHnzZe9r0QIGDkwkliuOFw/nXHlk9xrv0QNefx1uuw2qqjAJqqpg0CDo3TvppK4AXjyccw1v0SL4y1/gkkvCXOP33w+rrBIKxeTJjB01CiZP9sKRIl48nHMNK9Nr/I47Qq/xm2/2XuONgPfzcM41nOxe4zfeCMcfn3QiVyJePJxzDWP6dNhzT/jgA7jvPjjwwKQTuRLy4uGcK73334fddw9HHs8841PGNkJePJxzpfXaa9C1a2jXGDMGfvvbpBO5BuAN5s650nnuuXCU0aoVvPiiF45GzIuHc6407r0XunWDzTcPhWPzzZNO5BqQFw/n3PK79trQR2PHHWHsWFhvvaQTuQbmxcM5V39mcM45cOqpoS/Hs8/6XONNhDeYO+fqZ9GiMMf44MFhXg6fMrZJKfuRh6StJY2UVC3pM0kXSKrzHSepg6TnJH0dbyMkbVeOzM65HPPnQ8+eoXD8/e9wyy1eOJqYshYPSa2BEYAB+wEXAKcD59ex38Zxv2bAn+KtGTBcUlVDZnbO5Zg9O/TheOIJuOEGOP98kJJO5cqs3KetjgWaAz3MbC7hw391YICky+K6fLoBqwHdzewbAEkvAbOArsDNDR/dOcdnn8Eee4ROgEOHwkEHJZ3IJaTcp632AoblFImhhILSsZb9VgIWAd9lrZsX1/lXHufK4YMPYIcdwui3zzzjhaOJK3fx2BKYmL3CzKYC1fG+mjwUt7lS0rqS1gWuBmYDDzRQVudcxn/+Ey7Dra4OvcZ32SXpRC5hspqmgmyIJ5MWAmea2TU566cB95jZObXs+2vgSWDDuGoGsJeZvVXD9n2BvgBt2rRpP3To0HplnjdvHi1btqzXvklIU940ZYV05S1l1tbjx7PteefxQ+vWvH3ZZczfaKOSPG62pvralsPy5O3cufMEM+uQ904zK9sNWAicmmf9NODiWvZbH/gQeAzYM96eiPu1ret527dvb/U1evToeu+bhDTlTVNWs3TlLVnWe+81W2kls1/9yuyzz0rzmHk0yde2TJYnLzDeavhcLfdpq9lAvh5EreN9NTmT0O5xgJk9a2bPAj2BxcAZJU/pnIPrroNDDw3tHGPHwvrrJ53IVZByF4+J5LRtxMtwW5DTFpJjS+AdM1uYWWFmPwDvAJs1QE7nmi4z6NcPTjkFunf3XuMur3IXj2eAPSStlrWuFzAfGFvLflOAbSWtnFkh6WfAtsDkBsjpXNO0aBEcfTRcfHHoNf7AA2GucedylLt43AIsAB6WtGts1B4AXGVZl+9KmiRpcNZ+twMbAI9I6iZpb+BRQlvIoHKFd65Rmz8fDjgg9Bo/7zzvNe5qVdZOgmY2W9IuwA2EBu85hEtuB+TJtWLWfhMk7Qn0B/4ZV/8X2M1quNrKOVeEOXNg333hhRfg+uvhxBOTTuQqXNkHRjSzd4Fa56Q0s3Z51o0ERjZQLOears8+C3ONT5zovcZdwXxUXeeasg8+CMONzJoFTz8Nu+6adCKXEl48nGuqxo+HvfYKgxqOGQPt2yedyKWITwblXFM0fDh07gwtW4YpY71wuCJ58XCuqRk6NMw1vumm8NJL8POfJ53IpZAXD+eakuuvD73Gt9/ee4275eLFw7mmwAzOPRdOPhn22w+GDYM11kg6lUsxbzB3rrFbtAiOOw5uvz30Hr/pJmjm//Xd8vEjD+cas/nz4cADQ+E491y49VYvHK4k/F3kXGMyZAj060fHqVNhww1h1VVDX47rroOTTko6nWtEvHg411gMGRIGM6yuDnMzT5sW1p94ohcOV3J+2sq5xqJfvzBNbK4nnih/Ftfo1Vk8JI2StGX8/XBJazV8LOdc0aZOLW69c8uhkCOPPwJrxN/vxCdfcq4yrbNO/vVt25Y3h2sSCmnz+BQ4UNI8QMAm8fe84qi5zrlyuu++MLihFPp0ZLRoAQMHJpfLNVqFHHn8AzgZeAsw4F7CXBq5t//Fn865crrhBjjkENhxx3ApblUVJkFVFQwaBL17J53QNUJ1HnmY2W2SHgd+DowDTgD86MK5pJnB3/8OF10Ueo3/+9/QvDkcfTRjx4yhU6dOSSd0jVidxUPS4cBTZvaCpPOBx8zss4aP5pyr0aJFcPzxcNttcNRRYcpY7/znyqiQ01bZjeR/BzZquDjOuTp9/33oNX7bbeHy3Ntu88Lhyq6Qd9xsYIP4uwjtHs65JMyZE05RjRvnvcZdogopHiOAf0p6Py7fJem7mjY2s9+XJJlzblkzZoS5xt97L7RvHHxw0olcE1ZI8TgSOA7YEvgt8AnwZUOGcs7l+PDDMNf4zJnw1FOw225JJ3JNXCFXW1UDVwJI2hXoZ2Zv1fcJJW0NXA9sD8wBbgfON7PFBezbAzgb2BaoBv4D9DSzGo+EnEu9CRPCXONmMHo0/O53SSdyrriBEc1sk+V5MkmtCafB3gX2IzTEX0louD+3jn3/AtwAXAacCbQGuuCDO7rGbMQI6N4d1loLnnsOttgi6UTOAfX44JW0KeHDeydgTeBr4HngCjP7uI7djwWaAz3MbC4wXNLqwABJl8V1+Z5zbeBq4CQzuy3rrkeKze9catx/Pxx2GGy5JTz7LGywQd37OFcmRY2qK6k98CbQk3DK6J74syfwhqTf1vEQewHDcorEUEJB6VjLfgfFn3cXk9e51LrxxtAgvt124coqLxyuwhQ7JPsVwBtAOzM70szONrMjgU3i+ivq2H9LYGL2CjObSmi/2LKW/bYD3geOkjRN0kJJr0raocj8zlW2TK/xE0+EffYJp6p8rnFXgWRWeLeNeInuQWb2VJ779gbuM7NVa9l/IXCmmV2Ts34acI+ZnVPDfsOAHYC5wFnAV/FnB+DnZvZFnn36An0B2rRp037o0KEF/Y255s2bR8uWLeu1bxLSlDdNWaEMeRcvZotrrmGDJ59kxl578cHpp2Mrrlivh/LXtuGkKSssX97OnTtPMLMOee80s4JvwCzg8BruOxz4qo79FwKn5lk/Dbi4lv2eI3RO3DNr3eqEDowX1pW7ffv2Vl+jR4+u975JSFPeNGU1a+C88+ebde9uBmZnn222ZMlyPZy/tg0nTVnNli8vMN5q+FwttsH8KeASSR+b2QuZlZJ2Ioy+W9eUZbOBVnnWt4731bafAWMyK8xsrqQJwNaFRXeuQn3zTeg1PnYsXHMNnHJK0omcq1OxxeP/gMeAsZJmAjOBdePtZeD0OvafSE7bhqSNgRbktIXkeI8wNIpy1gtYUmh45yrOjBmhD8c774Q5yA89NOlEzhWkqAZzM/vKzHYCugE3AS/Gn3uZ2U5m9lUdD/EMsIek1bLW9QLmA2Nr2e/J+LNzZoWkVkB7wjwjzqXPpElhDo5Jk0KvcS8cLkWKOvKQtKKZLTazZ4Fn6/F8txAmlnpY0qXApsAA4CrLunxX0iRgrJkdBWBm4yU9BgyW9DdC28tZhDaUG+uRw7lkvf56GKdqyRIYNQp+70PCuXQp9lLd6ZIuk7RVfZ7MzGYDuwArEtpHzid0/uufs2mzuE22w4BHgauABwmFo0t8TOfSY+RI6NgxTNz04oteOFwqFdvmcQvhqqrTJY0HBgNDrYae4flYmOO8Sx3btMuzbh5hgMbjignsXEV54IHQa3yLLUKv8Q03TDqRc/VSbJvHADPbFNiN0GnvKmCGpCFx0ETnXE1uugl69QoDG44b54XDpVqxp60AMLNRZnY4sB5wEvD/gGGSJksaIMnHUnAuwwz694cTToC994bhw6F166RTObdc6lU8snQAdiZcfjubMEDiX4BJkg5bzsd2Lv0WL4bjjoMLLoAjj4SHHw5tHc6lXNHFQ1KVpP6SPgJGAusTJozawMz+BFQBtwKXlzSpc2nz/fdw0EFw663wt7/B7bf7XOOu0Sj2Ut3RwB+B6cCdwJ1mNiV7GzNbLOlewLvJuqbrm29g//1hzBi4+mo49dSEAzlXWsV+DZoJdAWGx3FPavImYaRd55qezz8PfTjeeQf+9S/o3TvpRM6VXLEzCfYqcLuFwJQ6N3SusZk0Kcw1/vnn8MQToYg41wjV6wSspI2ALYBVcu8zs6eXN5RzqfT662GcqsWLQ6/x7bZLOpFzDabYNo/VgPuB3TOr4s/sU1j1m4DAuTQbNSq0cbRuDcOGhaljnWvEir3a6h9AW0KjuYDuQCdCT/NPgD+UMpxzqfDgg+GIo23bMNyIFw7XBBRbPLoCA4FX4/JnZjbOzPoShmo/s5ThnKt4N98cLsft0CH0Gt9oo6QTOVcWxRaPNsCnZrYY+A5YM+u+p1l6Osu5xs0MBgyA44+Hbt1Cr/E116xzN+cai2KLx6fA2vH3D4G9s+7bDvi+FKGcq2iLF4eicf750KcPPPIItGiRdCrnyqrYq62GA7sCjxCGUr9bUntgAWGYkitLG8+5CvP992FU3Icegr/+Ff7xD1DuBJfONX4FFQ9JzQntHTOBjyW1MbN/SpoHHAA0B04kDEviXOMyZAj060fHqVNh5ZVhwQK46io47bSkkzmXmDqLh6RNgRFAu6zVcyUdZGaPEI5CnGuchgyBvn2hujpcl75gQSgg666bdDLnElVIm8dlwBLC5bktgG2AN/CjDNcU9OsH1dXLrvvhh7DeuSaskOKxPXCumb1oZt+b2XvAMUBbSes3bDznEjZ1anHrnWsiCike6wMf56z7iNBJcL2SJ3KuUoweXfN9bduWL4dzFajQS3VrG0HXucbnwQfDoIbrr//TyZtatICBA5PJ5VyFKLR4DJM0M3MDZsT1I7PXx/ucS7dbblnaa/y//4XbboOqKkyCqioYNMiHWXdNXiGX6p7f4CmcqwRmYbrYAQNCr/H77w9HGb17Q+/ejB0zhk6dOiWd0rmKUGfxMLOSFg9JWwPXExri5wC3A+fHIU8K2X8F4DWgPbCPmT1ZynyuiVq8GE46KYxVdcQR4WhjpZWSTuVcxSrrhMqSWhP6jLwL7AdsRuiVvgJwboEP8xfAR59zpbNgQeg1/uCDcNZZcMkl3mvcuTqUtXgAxxJ6o/cws7nAcEmrAwMkXRbX1SgWn4HA3whHLM4tn7lzwzwco0fDFVfA6acnnci5VCh2YMTltRcwLKdIDCUUlI4F7H8h8CIwsgGyuabmiy+gUyd4/nm45x4vHM4VodxHHlsCo7JXmNlUSdXxvidq2lHSL4EjgV82aELXNHz8Mey+O8yYAY8/HiZzcs4VTGbl68IhaSFwppldk7N+GnCPmZ1Ty75jgVfN7CxJ7QgzF9bYYC6pL9AXoE2bNu2HDh1ar8zz5s2jZcuW9do3CWnKm1TWlpMm8cuzzkKLF/Pff/yDuVtvXdB+/to2nDTlTVNWWL68nTt3nmBmHfLeaWZluwELgVPzrJ8GXFzLfgcDnwOrx+V2hI6LexfyvO3bt7f6Gj16dL33TUKa8iaSdfRos9VWM9t4Y7N33y1y19ENEqkhpCmrWbrypimr2fLlBcZbDZ+r5W7zmA20yrO+dbzvJyStBFwOXAqsIGkNYPV496qSVmuAnK4xevhh2GOPMFXsiy/CVlslnci51Cp38ZhIaNv4kaSNCaP1Tqxhn1UJl+ZeRSgws4G34n1DCSP8Ole7W2+FAw+E9u3hhRdg442TTuRcqpW7wfwZ4ExJq5nZt3FdL2A+MLaGfeYBnXPWrQf8GziHnAZ455ZhBhdeCP37Q9eu8MADPmWscyVQ7uJxC3Ay8LCkS4FNgQHAVZZ1+a6kScBYMzvKzBYBY7IfJDaYA/zXzF4tQ26XRosXwymnwI03wuGHw+23e69x50qkrKetzGw2sAuwIuGy3PMJc6H3z9m0WdzGufpZsAAOOSQUjjPPhLvu8sLhXAmV+8gDM3sX6FLHNu3quH8yYT4R535q7lzo3h1GjYLLL4czzkg6kXONTtmLh3MN6osvQtvGW2/B3XeH01XOuZLz4uEaj48/DpfiTp8eeo137Zp0IucaLS8ernF4660w89+CBTByJGy/fdKJnGvUyt3Pw7nSGzsWdt4ZmjULfTi8cDjX4Lx4uHTL9BrfcEN46SUocJwq59zy8eLh0mvQoNBr/De/CcOqe69x58rGi4dLn0yv8WOOCUcdI0bAWmslncq5JsUbzF26LFkCJ58cOv/96U8weLB3/nMuAX7k4dIju9f4GWd4r3HnEuRHHi4dvv029BofOdJ7jTtXAbx4uMo3c2bo8Pfmm95r3LkK4cXDVbZPPglzjU+fDo89Bt26JZ3IOYcXD1fJvNe4cxXLG8xdZRo3znuNO1fBvHi4yvPoo+FU1QYbeK9x5yqUFw9XWW6/HXr2hF//2ucad66CefFwlcEMBg6Eo48OvcZHjvRe485VMG8wd8lbsiTMNX7DDXDYYXDHHd75z7kK58XDld+QIdCvHx2nTg2npTbYAF55BU4/HS67DFbwA2LnKp0XD1deQ4ZA375QXR0moZ86NdwOPhiuuCLpdM65AvlXPFde/fpBdfVP17/8cvmzOOfqrezFQ9LWkkZKqpb0maQLJK1Yxz6/k3SnpElxv/cl9Ze0SrlyuxKZOrW49c65ilTW01aSWgMjgHeB/YDNgCsJRezcWnbtFbe9FPgQ+CVwYfzZswEju1Jbbz2YMeOn69u2LX8W51y9lbvN41igOdDDzOYCwyWtDgyQdFlcl88lZjYra3mMpO+BWyVVmdmUBs7tSmHcOJgzB6RwaW5GixbhMl3nXGqU+7TVXsCwnCIxlFBQOta0U07hyHgj/tygdPFcg8n0Gq+qgmuugaoqTArLgwZB795JJ3TOFaHcxWNLYGL2CjObClTH+4qxPbAE+Kg00VyDye01fvLJMHkyY0eNgsmTvXA4l0Ky7NMHDf1k0kLgTDO7Jmf9NOAeMzunwMdZD3gbeNrM+tSwTV+gL0CbNm3aDx06tF6Z582bR8uWLeu1bxIqKq8ZbYcMYdPBg/lqu+14p39/ljRv/uPdFZW1AGnKm6askK68acoKy5e3c+fOE8ysQ947zaxsN2AhcGqe9dOAiwt8jJWBccDHQOtC9mnfvr3V1+jRo+u9bxIqJu/ixWYnnWQGZocdZvbDDz/ZpGKyFihNedOU1SxdedOU1Wz58gLjrYbP1XI3mM8GWuVZ3zreVytJAu4BtgF2NLM693EJWLAAjjgC7rvPe40710iVu3hMJKdtQ9LGQAty2kJqcA3hEt/dzKyQ7V25ffst9OgBI0aEonHmmUkncs41gHIXj2eAMyWtZmbfxnW9gPnA2Np2lHQ2cCJwkJm90LAxXb3MnBmmiX3jDbjrrnD04ZxrlMp9LuEWYAHwsKRdY6P2AOAqy7p8N/YkH5y1fChwMeGU1XRJf8i6rVPeP8Hl9cknsNNO8M47Ya5xLxzONWplPfIws9mSdgFuAJ4A5gBXEwpIbq7sIUt2jz/7xFu2PwN3lTSoK87bb4c5OBYsCKerdtgh6UTOuQZW9lF1zexdoEsd27TLWe7DT4uGqwTjxsG++8Jqq4UJnHzKWOeaBL8ExtVfptf4+uv7XOPONTFePFz9+FzjzjVpXjxccXyuceccPpOgK0b2XON/+hMMHuxzjTvXRPmRhyvMggVw6KGhcJx+eujH4YXDuSbLjzxc3bJ7jV9+OZxxRtKJnHMJ8+LhajdzJnTtCm++CXffDYcfnnQi51wF8OLhavbJJ+FS3OnTQ6/xbt2STuScqxBePFx+b70Fe+4Z2jpGjoTtt086kXOugniDufupsWNh552hWbPQh8MLh3MuhxcPt6xHHgn9Nzbc0HuNO+dq5MXDLXXbbXDAAfCb38Dzz3uvcedcjbx4uNBr/KKLoG/fcNQxYoT3GnfO1cobzJu6JUvg5JPhxhu917hzrmB+5NGULVgAhxwSCscZZ3ivcedcwfzIo6n69lvo3j1chuu9xp1zRfLi0RTNnAl77RX6cnivcedcPXjxaGqye40//ngYesQ554rkxaMp8V7jzrkS8QbzpsJ7jTvnSsiLR1PgvcadcyXmxaOxGzTIe40750qu7MVD0taSRkqqlvSZpAskrVjAfq0k3SlptqRvJA2R5N2ga2IGF14IxxwT2jm817hzroTK2mAuqTUwAngX2A/YDLiSUMTOrWP3+4EtgL8AS4BLgUeBPzZQ3PRavDjMNX7jjeEy3Ntv985/zrmSKveRx7FAc6CHmQ03s1uA84H/k7R6TTtJ2h7YHTjCzB4ys0eAw4CdJO3aIEmHDIF27ejYpQu0axeWK1l23tVXX9pr/M47vXA450qu3MVjL2CYmc3NWjeUUFA61rHfF2Y2LrPCzF4DPon3ldaQIWGQwClTkBlMmRKWK7WA5Oatrg4F49e/hhW8Wcs5V3rl7uexJTAqe4WZTZVUHe97opb9JuZZ/168r7T69QsfwNmqq6FPH7j44pI/3XL74ANYtGjZdQsXhr+jd+9kMjnnGrVyF4/WwJw862fH++qz36b5dpDUF+gL0KZNG8aMGVNwyI5Tp6I8623RIr5cZ52CH6dc1nn33fx5p05lbBF/d7nNmzevqH+XpKUpb5qyQrrypikrNFzeRtvD3MwGAYMAOnToYJ06dSp857Ztw6mqHKqqYt1KfNO0a5c/b9u2FPV3l9mYMWMqOl+uNOVNU1ZIV940ZYWGy1vuE+KzgVZ51reO95V6v/oZOBBatFh2XYsWYX0lSlte51zqlbt4TCSnjULSxkAL8rdp1LhfVFNbyPLp3Tt0rquqwiSoqgrLldp+kLa8zrnUK3fxeAbYQ9JqWet6AfOBsXXst56knTIrJHUgtHc80xBB6d0bJk9m7KhRMHly5X8Qpy2vcy7Vyl08bgEWAA9L2jU2ag8Arsq+fFfSJEmDM8tm9jLwHHCPpB6S9geGAC+Y2Yhy/gHOOefKXDzMbDawC7Ai4bLc84Grgf45mzaL22TrRTg6uQO4B5gAdG/IvM455/Ir+9VWZvYu0KWObdrlWTcH+HO8OeecS5B3P3bOOVc0Lx7OOeeKJjNLOkODk/Ql8NNedIVZG5hVwjgNLU1505QV0pU3TVkhXXnTlBWWL2+VmeUdVqNJFI/lIWm8mXVIOkeh0pQ3TVkhXXnTlBXSlTdNWaHh8vppK+ecc0Xz4uGcc65oXjzqNijpAEVKU940ZYV05U1TVkhX3jRlhQbK620ezjnniuZHHs4554rmxcM551zRvHg455wrmhcP55xzRfPi4ZxzrmiNdg7zxi7OwNgVEPCAmX0laSPgDGAzYDIwyMz+m1xKkPRX4OmkcxRKUnOgmZl9m7VuHeBEYGtgCfAmcJOZfZNISOcqgF+qG0kSYX6QbsBWwJqED4rPgVeAu8zsg+QSLiXp98BwYFVgEfA1sAfwNLAYeAfYFlgP2NXMnk8oKpKWAEaYLvheYKiZfZRUnrpIehr40MxOicvbE2arXEKYQ0ZAe+AHoIuZvZNg1t8Azc3spax1ewJns7TQvQUMyN6mUsT/c/sAvyW8R8YTvmhU9IeSpNUJY0V1MbMXks4DP2bqAqwMPGVm38UvPScQZlz9mPBl8rOSPWeF/zuVRXyRnyZ8KHxO+GDYkPCGfobw4v8/4EIzuzCpnBmShhOOGrsD3xEm1Nqf8OF2gJktlPQz4FFgFTPrnFDUTPG4FPgFsBsh9+uEQnK/mU1PKls+kmYBR5nZY3H5FcJrvH/maERSK+Bx4Hsz2yPBrK8AT5jZwLh8JHA7MBoYRSh0uwB/BHpm/qaEsr5EeF3fi8utCbODtgfmxc1aEr6o7ZF95JcEScfXcndz4HLgWuBDADO7qRy58pG0OTAS2Diu+gTYnfAFcw3gI8Ln13ygvZlNK8kTm1mTvwH/JrwJfpG1bgPgWeChuNyR8CY/sgLyfgXslbW8LuFb5u4523UDZiWcdQnw+/j7mkDf+EZfFG9j4ro1k35dY8ZqYOes5R9yX9es1/a7hLPOzc4GTAKuz7PdLcBblfI+iMuDCUfMe2at2xOYDVxdAe+DJYSj+CU13LLvW5xw1vsJR5ibx/9j/4yfZy8Bq8Vt1o7b3Fqq5/UG82Av4G+WdV7ewuHdscD+ktY3s7HAxcApCWXMZXl+zz2MrKjDSjP72swGmdkuwEbA6YTD7FuAGZKeSjRg8D8g+0jtC8J/yFxrEQpNkpbkLFcBD+bZ7kHCN89Ksi9wgZk9m1kRfx8I9Egs1VKPATOBo4AVzWyFzI3wfhDQKa7LnTK73HYCBprZJDP7GjiX0O55hcUjODObBVzDsu/t5eLFI1iB8E0i12LCm6RVXH4V2KJcoWoxHjhDUktJKwDnANOB4yStCCCpGXA84cOw4pjZ52Z2rZntAGxCmMd+g4RjAVwC/E3SkfE1HAhcLmk3SStL+llsV/gH4Rtfkp4HemctvwPkG3r7d4T3RyVZg/A+zjWB0FaXKDPrDhwBnAn8R9KO2Xcnk6pGrQmn2zMy/9a5cxh9TPjSVhJ+tVUwHLhI0ttm9jH8eE72OsI/SqahvCVQCVfY9CNknk049VNNaCx7APhQUqbBfAPCqYCKZmZTCB/al1RAloclnUT4lnY18D7hy8OzhA8NxU0fJ3ywJOkc4MX4BeJ6QkP53ZLWJJwOFOF9cSrwt4QyZuspKVPcZgP5Jhlam3A6LnFm9pykXxJev6ckPUu4mjHR9pg8ZhKOOjMWA7cSjpqzrUsJs3uDORAvcR1GOKqYQjjPvQmwADjEzJ6J211GmFmrV1JZM2LmvQlfAB4ysxmS1gPOIpyimALcbmavJxgTSf2B26yEV3mUg6S1gF7A7wnfhEX4wHsPeNLMJiQY70eSfg3cDGzHssUt8/tswumhaxMJGMULJ3LdZWZH5mx3K7C1mf2xPMkKE/9vXUY4pXYroaB0NrNxSeYCkPQo8HXua5lnu+uBLc1st5I8rxePIJ7uOQj4FbAKofHx3ngO0bmKJmkrQgHJLXQvmdnCJLMVQ9LRwEdmNirpLPnES7evJnxB62YVcAm0pDZACzP7pI7t/o9w4cTIkjyvF4/GR9IKZpbvm17FkLQKoVFvCTCpEj/gYpvHpmT1+TGzqcmmcq4yeIN5DknbSOop6S/x1lPSNknnyiWph6RHJT0taZ+4rpekycBCSVPit7hESTos9j/ILDeTdAnhMs23CQ36X0uqhHPyAEhqL+lxwvnh94AXgZeBTyRNl3SBpBaJhmxEFCWdIx9JzXP/rSX9On4utE8qV0VI8vrkSroBRxLaCfJd272YMNzHn5POGbMeFHO9QLiksBo4mtBWM5jQq/TfMfceCWd9Fzgua/nKmPc8YEfCZYYDCB2YzqmA13Z3QlvXeMKl2QMIHUV/iJlPJ1zV9CbQugLy7k3oN/NefC/snGeb7Ui+L8LuxD4HWev2J3QYXQQsjK95t6Rf05itFfBIzLUIuA1YEbg753PhBWDtpPMW+Df1LOX7IPE/qBJuwEnxTXIjoTfu2vGNsmL8fSfghvihckIF5P0PcEvWcu+Y7cqc7e4ERiSctRromLU8Ezglz3ZnAFMq4LWdANxdw3tkMuFofZX4oXdTwll3ix9gL8b354S4fCXxlHTcrhKKx2KW7STYPX4AvxT/7U8nHN0tIk+nzATyXkcYguQk4PD4heEh4NNYCNch9A+bDtycdN4C/6aSFg9v8wAkfUz4ML6sju3OAo41s03Lk6zGHHOBHmY2Ii63IjSQ7mpZDY3xdNatZpZY/wlJM4ATzeyhuLyAcDQ0Jme73YDHzax5+VMuk2M+sK+ZDc9Z35rQs38bM3tP0uHApWa2fhI5Y6YXCONw/Tlr3ZGED77hhCsFv5e0HaHhPLHObPFqqz+Y2Wtx+XVgupntk7Pd08CqZtYxgZjZOSYTOt7dFpd/QyjOfzazu7O2O5pwxLxJIkFDhjsK3LSK0LGxJO8Db/MI1gdeK2C716iADkyEyzCz3wCZsYHm5Gw3j9AZK0mPEzo0rhyXRwCH5NnuEMK3u6TNJFxxl+tXhNc9089nCks7jyZlW+Bf2SvM7A7CUDp/AEbFPh+VaFvCJa+5BhEGSkzaOizt3wVxDCvCOFHZJpG/v0o5HUE4GvpFHbeqmh6gPryTYPAWcLSkcVbDVUqxQe9oQiNv0qYQRk0dBmBmi+MlhO/lbLcZMKPM2XKdTegJ/T9JtwNPAJdK2palHdk6A78hjLCatEHAhZJaEgrdD4Qe2v2A0ba0v8qmQNJXXn1PGFl5GWY2IfaIHkY4LTSgzLlqkn2a4xvyd1j7jsr4UvsJoQiPjct/JJxm24HQzpGxI8m/Dz4EXjOzw2vbSNIBwH2lelIvHsHphB7E70p6mDB8+Jx4XytgS8I52o2ojB7bD5MzzICZvZpnu4NZ9o1edmb2taQ/ED58/4/QyxVg+3j7gXCK5Y9m9p9kUi5lZgPjKZa/AX/PrCZcgHBq1qYLCQ3qSXqbcN798dw7zOzjWECeBu4qc66aDJO0KP7eivCFYWzONluS/BceCOOtXSvpF4RCdxDhi9B58YvFW4QjpNOApEfafoVQ1OqS3Yl0uXmbRyRpM0Lv7D1ZOrRxxqeEK24utwqeiyKXpLbAHDOriOEeACS1Y9mObB9ZZfbxWIlw5LYK8HElvYYZko4hDFHyG6uhM6ukVQlXDe1qYVC/RMSRBnJ9aGb35mw3Jq6vhMvMTyacTl2JMFrDLZIOIbQpZQbGHAT8Ncn3cLxkeEczu66O7dYmtNnlFuz6Pa8Xj5+K13WvERfnmFnSo6c65ypEPIW9tpl9mXSWJHnxaGTiIfXrQO9KOA2kFE7rqpRM8etckrx4ZIkfGhsQTqXMynP/2kBXM7un7OGWzdG1lrtXJTSK/Y04HLuZPV2OXPkoRdO6Qrqm+C1UHPfqQDO7IOEcZZ8qtZTiEUf2tLkTCH9H4h+iCqMV9yT8f7rLzCZK+hVwPku/8NxgZsNK9qRJd1yphBvwM8Jw5ovjbSGhp3arnO0S72wVc6RplrNZwH5Zy68QekSvlrWuFaHhdFgFvLbDCdO4rkE4130DMI3Qe3ulrPfLM4SrrxJ//xbwN5W0c1g9M2xOuEow8778iPCh9jGhQP+HMBT7F8BGFfCavQRslbXcOmZcEnPOZWknx9WSyhmz7UH48vV5fF3nEq5gnB3z3Rj/3y0mTKdcmudN+h+pEm6Eq2rmEC7F7QCcHN/EHwI/z9quUorHeMIVKX8mXLudfftlfFMflFmXcNbUTOsac6Rpit+2Bd6OTfp9S0JTpS5H3tRMm0sYYeABwoyHEC6imA0Mztnun8ArJXvepP+RKuFGuDT3xJx16wHjgC+B7eO6SikeIsz7PZMwZMImWfe1im/8n4xxlFDW14D+WcufAgfn2e5w4MsKyPtVzgfEOvH13C1nu64VUDwyR5l13SrhCPQz4KCs5aqYq0fOdn8GPqiA90Fu8fgSODXPdokPq0O4lHjXrOXWMX+XnO12J1wAVJLn9X4ewcbkdP4zs88l7UKo1iMk9aYyrj/HwjthkKT7gYuA/8aJXi5KNllelwBDJH0K3MPSaV2/IpyqynQSrIRpXWHpFL8vEI6asqf4HWWhQ2alTPH7LTAKuL2O7XYiXIaepESmSi2hNajcaXPns2xn0czvuUP9tCB0LC0JLx7BZ8DPCUcaP7Jw7fbBkq4lHBYm2lCey8zmACdKGkS49vxD4FIqaI5lS9e0rpCuKX5fI7TLPVXbRnHulKQlMlXqckrLtLkvAn+X9GHMcgVhNOu/xlEzvo3j351FKHYl4Vdb8ePAYpuaWadatjmb8K3ZLMEB5moj6WDCVJkbEQZAS3yKzAylZFpXSNUUv+cBfc0st1Nr7nY7A+ebWefyJMub4VESmCq1vpSiaXMlbU4YSifzPphMOJp/kNBjfwrQjvBlqLOZvVmS5/Xi8eNlbr2Af1gt085KOpRw7vvPNW2TtHhKZVVgnpktTjqPc5DcVKkNTRUybW7s37Uj4QrBkWY2P3Z2/gtLv/Dca2bTSvacXjycc84VqxJGr3QNRNJtkgYnnaMQacoK6cvrXKl5g3kRJN0GrGBmRyWdpUCdSc8XhDRlhRTllTSCcJZhl6Sz1CVNWSFdeUud1YtHcVLzgQFgZpsnnaFQacoKqcsr0vO+TVNWSFfekmb1No9GLF6iua6ZJT1ZTZ3SlBXSl9e5UktLxawIklaJc2SkRTfCjGhpkKaskKK8klZKy/s2TVkhXXlLndWLR3FS84HhmgZJJ0j6SNJ8SW9J+lOezX5LBbxv05QV0pU3iaze5pFCkgq9pjxfj9iySlNWSFfe2Cn0esIUuW8QpiK9S9J+wGFmVrKhKJZXmrJCuvImldWLB+n6wIh2Jgzz8W4d21XCsBRpygrpynsGcIWZ/ThuVRyPbQgwWtLeZvZVYumWlaaskK68iWT1BnNA0iIK+8DYENgu6eFJJL0FTDSzXnVsdwBwX5J505Q15khNXknfAvuY2Zic9e0I842sSBh/ax3gJc9auDTlTSqrt3kE7wD/M7MDa7sBVyUdNHoF+EMB22UPPJiUNGWFdOX9hjAw3zLMbDLh1MUs4GXgd+WNlVeaskK68iaS1Y88+HFwsz3NrKqO7XoS5rROtOhK2gzYxswer2O75oTLSXOHvS6bNGWNOVKTV9JjwLdmdlgN9zcnDI63FwkP6JmmrDFPavImldWLB+n6wHAuQ9KBwGnA3jUN6ClpReBmwoCem5QzX06O1GSNWVKTN6msXjycc84Vzds8nHPOFc2Lh3POuaJ58XBNjqQ+kiZI+lbSbElvSGqQK+kkbSFpgKQ1Cth2gCTLun0m6aHYJlfXvndJyjfHtnMNwouHa1IUphO+HRgG9AAOBx4D9m2gp9wC6A+sUeD23wDbx9sZwK+BkZJWrWO/C4E+9UroXD14D3PX1JwI3Gpm52Ste0LS+UkFyrHIzF6Jv78iaSrwPNAVeCB3Y0nNzWy+mX1UzpDO+ZGHa2rWAD7PXWlZlx1KahdPGx0q6Z/x9NZMSf1z95PURdKrkr6X9IWkmxTmk0ZSJ+CJuOkn8TEnF5l3QvzZLj7mZElXSjpP0jRgblz/k9NWkqok/VvSLEnVkt6WdGjW/atIukzSp5IWxAH1uhaZzzVRfuThmprXgZPiN/on6xjz53LgSeAAwphX/SXNMrMbASRtAzwLDAd6AhsDlwCbEoaDeJ047hDhFNkMYEGRedvFn9kF71DCqAjHU8P/YUnrEnoVV8cMnwLbxowZDwK/J5xW+wg4CHhcUgcze7PInK6J8eLhmpoTgEeBuwCT9B7wEGFgubk5275jZsfE34fFD+RzJN1sZkuA84ApwL5mthhA0tfAfZK2N7OXJb0f938jDhdRJ0mZ/5ebAjcB3wIjcjbbu47RUk8DWgHtzWxGXDcy6zl2IUwx0MnMxsbVz0naAugHHFhIVtd0+Wkr16SY2dvAVoQG8psI41OdB4zPnG7K8kjO8sPABsBGcfn3wCOZwhE9BCwCdqpnxLWAhfH2PqGA9MoqAAAjCxhmuwvwbM5+2XYlHM28KKlZ5kYoMB3qmd01IX7k4ZocM1tAaIt4AkDSUYQrsI4Crs3adGbOrpnl9YGp8ecXOY+9WNJXwJr1jPcN4YPdCB/un9lPh4H44id7/dRawH9quX9tYD1Ckcq1OM8655bhxcM1eWY2WNJlwJY5d61bw/KMrJ/LbBPHEFoLyDvGUAEWmVld/TUKGVPoK0Jxq8nXwHRg/wJzObcMP23lmpTYbpG7bh1C+0DuN/ruOcuZRu9pcflVoHssGNnbNANeiMs/xJ/lnjxqJLCHpDa13L8eMM/MxufeyhfTpZUfebim5r9xCOvnCKehqghXI1UDd+dsu00crv8hwtVWRwGnxMZygIsI034+KulmQlvIpcAwM3s5bpNpMD9G0lCg2sz+2zB/2jKuJnSAfF7SQMLVVlsBq5rZZYQrxIYBwyVdSrh6a3VCp8RVzOzsMmR0KebFwzU1FwD7AdcR2iU+B14iNEp/krPtWcDehOLxPaEX9w2ZO83sHUl7ARcTGtPnEuaRPitrmymSzgBOBk4iHLW0a4g/LJuZfSlpR+Ay4BrgZ8CHwD/i/SapB3AOcCrQlnAq603CfNjO1cqHZHcuR5y+8xPC1J5PJhzHuYrkbR7OOeeK5sXDOedc0fy0lXPOuaL5kYdzzrmiefFwzjlXNC8ezjnniubFwznnXNG8eDjnnCva/weH2pgYE2jh4QAAAABJRU5ErkJggg==\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.minimum(np.maximum(0, x - strike_price_1), strike_price_2 - strike_price_1)\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": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact expected value:\t0.5695\n",
      "exact delta value:   \t0.9291\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(\n",
    "    uncertainty_model.probabilities[np.logical_and(x >= strike_price_1, x <= strike_price_2)]\n",
    ")\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": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "problem = EstimationProblem(\n",
    "    state_preparation=bull_spread,\n",
    "    objective_qubits=[num_uncertainty_qubits],\n",
    "    post_processing=bull_spread_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": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "result = ae.estimate(problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:    \t0.5695\n",
      "Estimated value:\t0.5686\n",
      "Confidence interval: \t[0.5610, 0.5763]\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": [
    "### Evaluate Delta\n",
    "\n",
    "The Delta is a bit simpler to evaluate than the expected payoff.\n",
    "Similarly to the expected payoff, we use comparator circuits and ancilla qubits to identify the cases where $K_1 \\leq S_T \\leq K_2$.\n",
    "However, since we are only interested in the probability of this condition being true, we can directly use an ancilla qubit as the objective qubit in amplitude estimation without any further approximation."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "# setup piecewise linear objective fcuntion\n",
    "breakpoints = [low, strike_price_1, strike_price_2]\n",
    "slopes = [0, 0, 0]\n",
    "offsets = [0, 1, 0]\n",
    "f_min = 0\n",
    "f_max = 1\n",
    "\n",
    "bull_spread_delta_objective = LinearAmplitudeFunction(\n",
    "    num_uncertainty_qubits,\n",
    "    slopes,\n",
    "    offsets,\n",
    "    domain=(low, high),\n",
    "    image=(f_min, f_max),\n",
    "    breakpoints=breakpoints,\n",
    ")  # no approximation necessary, hence no rescaling factor\n",
    "\n",
    "# construct the A operator by stacking the uncertainty model and payoff function together\n",
    "bull_spread_delta = bull_spread_delta_objective.compose(uncertainty_model, front=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "problem = EstimationProblem(\n",
    "    state_preparation=bull_spread_delta, objective_qubits=[num_uncertainty_qubits]\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": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "result_delta = ae_delta.estimate(problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact delta:    \t0.9291\n",
      "Estimated value:\t0.9290\n",
      "Confidence interval: \t[0.9269, 0.9312]\n"
     ]
    }
   ],
   "source": [
    "conf_int = np.array(result_delta.confidence_interval)\n",
    "print(\"Exact delta:    \\t%.4f\" % exact_delta)\n",
    "print(\"Estimated value:\\t%.4f\" % result_delta.estimation)\n",
    "print(\"Confidence interval: \\t[%.4f, %.4f]\" % tuple(conf_int))"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2019-08-22T01:55:52.763931Z",
     "start_time": "2019-08-22T01:55:52.753702Z"
    }
   },
   "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_ibm_provider</code></td><td>0.6.1</td></tr><tr><td><code>qiskit_finance</code></td><td>0.4.0</td></tr><tr><td><code>qiskit_algorithms</code></td><td>0.2.0</td></tr><tr><td><code>qiskit_aer</code></td><td>0.12.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:15:46 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
}