{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Pricing European Put Options"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "<br>\n",
    "Suppose a <a href=\"http://www.theoptionsguide.com/put-option.aspx\">European put 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\\{K - S_T, 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\\{K - S_T, 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 \\leq 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://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": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:37:27.742240Z",
     "start_time": "2020-07-13T23:37:25.358136Z"
    }
   },
   "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": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:37:28.704565Z",
     "start_time": "2020-07-13T23:37:28.692784Z"
    }
   },
   "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": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:37:29.373028Z",
     "start_time": "2020-07-13T23:37:28.931666Z"
    }
   },
   "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 decreases linearly with an increasing spot price at maturity $S_T$ until it reaches zero for a spot price equal to the strike price $K$, it stays constant to zero for larger spot prices.\n",
    "The implementation uses a comparator, that flips an ancilla qubit from $\\big|0\\rangle$ to $\\big|1\\rangle$ if $S_T \\leq 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 scaling 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://arxiv.org/abs/1806.06893)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:37:29.827435Z",
     "start_time": "2020-07-13T23:37:29.821279Z"
    }
   },
   "outputs": [],
   "source": [
    "# set the strike price (should be within the low and the high value of the uncertainty)\n",
    "strike_price = 2.126\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]\n",
    "slopes = [-1, 0]\n",
    "offsets = [strike_price - low, 0]\n",
    "f_min = 0\n",
    "f_max = strike_price - low\n",
    "european_put_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",
    "european_put = european_put_objective.compose(uncertainty_model, front=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:37:31.421735Z",
     "start_time": "2020-07-13T23:37:31.147717Z"
    },
    "tags": [
     "nbsphinx-thumbnail"
    ]
   },
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAY8AAAE+CAYAAAB1DJw3AAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjUuMSwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy/YYfK9AAAACXBIWXMAAAsTAAALEwEAmpwYAAA6R0lEQVR4nO3dd5hU9dnG8e8NqDQFbBQpK0mMGjWvYo0FwUqxAApGLNhQYy+xEQMWLFijRmwYSzBYUKQpAiLGLmiMCWpsgDQVxAJroTzvH7+zMgyzZXZn58zZeT7XNdfunDln995hmGfO+TWZGc4551w26sUdwDnnXPJ48XDOOZc1Lx7OOeey5sXDOedc1rx4OOecy5oXD+ecc1nz4uEKmqQhkizltkDSaEm/iClPE0mjJC2J8gyItp8i6VNJKyW9UM6x+6b9LWW3lfn8G6IsW0XPbfO07QOiTE3zncklS4O4AzhXBd8AB0ffdwSuAqZK+o2ZLc9zltOBQ4DjgPnAx5JaAcOBO4DHgaWV/Iz+wCcp9+MYbLUVMBh4APg6ZfsEYA+gNP+RXJJ48XBJsNLMXou+f03SXOCfQHfCm3U+bQ18YGajyzZI2guoD9xvZv+uws/4t5n9p7YC1oSZfQl8GXcOV/j8spVLopnR15LoMtIdkj6QVBpdOvqrpI3Kdpb0WKZLSdFlm88lrRfd31TSg9ElqVJJL0jaOWX/2cBJwI4pl5yGEAoZwDupl7KyFf2+J9K2lV3q2i66XxLd7yvpbknfSJon6QpJ9dKO3UHSOElfS1om6Q1JB0jaFxgX7fZp9PNmR8esc9mqsuel7LmRdKOk86I8S6PLe82r81y4wufFwyVRSfR1EdCY8Kl/ENANuBzoytpnJCOAfSRtWbZBkoDjgb+b2Ypo8xjgIOBCoB/h/8c0Sb+MHu8FTATeJ1za2QO4Dzgjerx/tG1CJfnrS2qQcqvO/8NhwDLgCODvwJ+j78v+vq2Bl4HWwGlR9qeAdsBb0d8I0DvK3KuC3zWGip+XMn2B/YCBwMVAT+CaavxtLgnMzG9+K9gbMARYTLjE2oBwrX4a8C3QOsP+DYA9Ce0I7aNt9YC5wBUp+3WN9tkuun9wdL9zyj5NCJdw7k7Z9gAwI+137pv6syr4W8r2S79dHT3+AvBERT+bUDgNeChtv38Bo1Lu/wOYBzQqJ0vP6OeUpG0fEG1vmuXzMhv4GGiQsu1WYFHcryG/1c7N2zxcEmwCrEi5PxfoZ2YLASQdC5wP/IrwxlZmK2Cuma2W9DfgOElDLLyzDSAUgbK2h12BL8xsetnBZrZc0nhgrxz/PUcR3mjLLKjGz3gu7f4soH3K/a6Es6rvq/GzU2XzvEwzs9SeY7OAzSWtZ2vO7lwd4cXDJcE3wP6ET8CLgAVRAUBSL+AhQm+ny4CvCJdqngIapvyMvxEuaXWR9CbQhzWXboiO+SLD7/4c2DiXfwzwX6t5g/nXafd/Yu2/dxNgYQ1/B2T3vGTKJGAD1i7+rg7w4uGSYKWZzSjnsSOB183sD2UbJHVO38nMZkuaQjjj2JJwKesfKbssBDbP8PNbEgpSPvwArJ+2rUU1f9YSwht/TRXC8+IKkDeYu6RrBPyYtq1/OfuOIJxx/AEYY2Zfpzz2OuESyz5lGyQ1BnoAL+UsbcXmEboCpzqwmj9rKtBXUsNyHv8p+lre42UK4XlxBciLh0u6yYSeVIMk7S/pZkKPn0zGED7d70S4jPUzM5sEvAI8Kul4ST0JPasaATfUVvg0TwG/knRL9LcMZc3gyGxdATQDXpTUL/p5f5R0YvT4B9HXUyXtJmn7TD+kQJ4XV4C8eLikuxu4CTgHeBLoABydaUcz+xF4BvgMmJJhl8MJxehWQldfAV3N7KNchy4n3wRCu80RhELSgfB3VednfUBo0F5M6E78VPRz50SPzyG0+fQmdOkdl/knATE/L64wKWp3dK7Ok9SA8OZ5v5ldHnce55LMG8xdnSdpfeC3hDOSTQhnK865GvDi4YpBG+ANQpfTU81sXsx5nEs8v2zlnHMua95g7pxzLmtFcdlq0003tZKSkmodu3z5cpo0aVL5jgUiSXmTlBWSlTdJWSFZeZOUFWqWd+bMmYvNbLOMD8Y9uVY+bp06dbLqmjZtWrWPjUOS8iYpq1my8iYpq1my8iYpq1nN8pI2CWjqzS9bOeecy5oXD+ecc1nz4uGccy5rXjycc85lzYuHc865rHnxKM/IkVBSQueuXaGkJNx3zjkHFMk4j6yNHAkDB0JpKQKYMyfcB+hf3lIRzjlXPPzMI5NBg6C0dO1tpaVhu3POOS8eGc2dm91255wrMl48MmnfPrvtzjlXZLx4ZDJ0KDRuvO72I47IfxbnnCtAXjwy6d8f7rkHOnTAJGjXLvS4Gj4cZsyIO51zzsXOi0d5+veH2bOZ/vzzoa3j1Vdh882hRw/45JO40znnXKy8eFRVq1bwzDOwYgV06wZLlsSdyDnnYuPFIxtbbw1jx4ZxH4ceCt9/H3ci55yLhRePbO21F/z97+Ey1rHHwqpVcSdyzrm88+JRHUccATfdBKNHw4UXxp3GOefyzqcnqa7zzguXr269NYz/OO+8uBM551zeePGoiZtugnnz4IILQndeHwfinCsSXjxqon59ePhhWLgQjjkm9Mjaa6+4UznnXK3zNo+aatQo9MDq0AEOOwzefz/uRM45V+u8eOTCJpuEMSANGoQxIIsWxZ3IOedqlRePXOnYESZMgC++gJ49YdmyuBM551yt8eKRSzvvDI8+Cm+/Df36wcqVcSdyzrla4cUj13r2hDvvhIkT4Q9/ALO4EznnXM55b6vacOqpYTLFa64JDem+AqFzro7x4lFbrr46FJA//SmMATnuuLgTOedcznjxqC0SjBgRxoCcdBK0aQP77x93Kuecy4m8t3lI2lbSVEmlkhZIulJS/Soct7Ok5yR9Fd2mSNotH5mrbf31w/xX22wDvXvDO+/Encg553Iir8VDUgtgCmDAYcCVwAXAFZUc1y46rgFwbHRrAEyW1KE2M9dYs2ah8XyjjaB7d/jss7gTOedcjeX7zOM0oBHQ28wmm9ldhMJxvqSNKjiuB7Ah0MvMJpjZBKAX0AToXtuha6xt2zCIcNmyUEC+/jruRM45VyP5Lh7dgElm9m3KtlGEgtK5guPWA1YCy1O2LYu2Kdcha8X228OTT8IHH4RLWD/9FHci55yrtnwXj62BtSZ/MrO5QGn0WHlGR/vcJGlzSZsDtwBLgcdrKWvu7bdfaESfNg1OPNHHgDjnEivfva1aAF9n2L40eiwjM1sgqQswHjg72rwQOMjMvsx1yFp17LGh3WPQoLAOyDXXxJ3IOeeyJsvjp19JK4A/mtmtadvnAQ+Z2WXlHNcaeBGYBQyPNp8B7Aj8Ljp7ST9mIDAQoGXLlp1GjRpVrczLli2jadOm1Tq2XGZsdfPNtBk/nv+ddx4LDj00Zz+6VvLWkiRlhWTlTVJWSFbeJGWFmuXt0qXLTDPbOeODZpa3G/AFMDjD9uWEolLecTcDs4H1UratD8wBbqvs93bq1Mmqa9q0adU+tkIrVpj16GFWr57Z2LE5+7G1lrcWJCmrWbLyJimrWbLyJimrWc3yAjOsnPfVfLd5vE9a20bUDbcxaW0habYG/mtmK8o2mNlPwH+BX9RCztrXoAGMGgU77QRHHQVvvhl3Iuecq7J8F49ngIMkbZiyrR/wPTC9guPmANtJWr9sg6QNgO0IZyTJ1LQpjB8PLVtCjx7w8cdxJ3LOuSrJd/G4C/gReFLS/lG7xBDgZkvpvivpI0kjUo67D2gDPCWph6SewBigNXBPvsLXipYtwxiQVavCQlKLF8edyDnnKpXX4mFmS4H9gPrAOMIAwVuAwWm7Noj2KTtuJnAwYaDgw8BDhEtdB5hZ8uf8+PWvw1K2c+fCoYfC99/Hncg55yqU94kRzWwW0LWSfUoybJsKTK2lWPHbc08YORKOPBL694fHH4f6lU755ZxzsfDFoApJnz5w883w1FNw/vk+iNA5V7B8SvZCc+654fLVLbeEhaTOPz/uRM45tw4vHoXoxhvDKPQLLgiTKvbtG3ci55xbixePQlSvHjz8cFhI6thjoXVr2HvvuFM559zPvM2jUDVsCE8/DVtuCYcdBu+9F3ci55z7mRePQrbJJmEMyPrrhzEgixbFncg55wAvHoVvyy3DKPQvvwyj0JctizuRc8558UiEnXcO4z7eeSc0nq9cGXci51yR8+KRFN27w/Dh4TLW6af7GBDnXKy8t1WSnHIKzJkDQ4eGhaQuvzzuRM65IuXFI2muuiqMAfnzn0MBOf74uBM554qQF4+kkeDee2HBAjj5ZGjTBg44IO5Uzrki420eSbT++vDEE7DNNmE+rHeSP7Gwcy5ZvHgkVbNmMHFi+Nq9e7iU5ZxzeeLFI8natg0FZNky2GMPaNeOzl27QklJmN7dOedqiRePpNt+ezjjDJg/H+bNQ2ahR9bAgV5AnHO1xotHXfDII+tuKy2FQYPyn8U5VxS8eNQFc+dmt90552rIi0dd0L59dtudc66GvHjUBUOHQuPG627v0yf/WZxzRcGLR13Qvz/ccw906IBJ0K5dmI13+HB444240znn6iAvHnVF//4wezbTn38+tHW8+iq0agU9e8LHH8edzjlXx3jxqKtatoRnn4XVq8NCUosXx53IOVeHePGoy7baCsaODaPPDzkkdN91zrkc8OJR1/3ud2Gw4Ouvh0tbq1bFncg5Vwd48SgGvXvDrbfCmDFw3nm+kJRzrsZ8SvZicfbZYdqSm2+GDh3gggviTuScSzAvHsXkhhtC+8eFF4ZJFfv1izuRcy6hvHgUk3r14KGHYOFCOO44aN0a9tkn7lTOuQTyNo9i07AhPP00dOwIhx0Gs2bFncg5l0B5Lx6StpU0VVKppAWSrpRUv4rH9pb0pqTvJS2R9KykJrWduc7ZeGN45hnYYIMwBmThwrgTOecSJq/FQ1ILYApgwGHAlcAFwBVVOPZk4BHgGaAbcDLwIX7prXpKSmDCBFiyBHr0gO++izuRcy5B8v3GexrQCOhtZt8CkyVtBAyRNCzatg5JmwK3AGeZ2b0pDz1V64nrsk6d4PHHwwDCvn3DgML11os7lXMuAfJ92aobMCmtSIwiFJTOFRzXN/r6YG0FK1rdusFdd4WpTE47zceAOOeqJN/FY2vg/dQNZjYXKI0eK89uwAfASZLmSVoh6XVJv6u9qEXk5JPh8svh/vvhqqviTuOcS4B8F48WwNcZti+NHitPK+DXwJ+Ai4FDgOXAs5Ja5jhjcbriCjj+eBg8GB54IO40zrkCJ8vjZQpJK4A/mtmtadvnAQ+Z2WXlHPcccADQzcyejbZtBMwB7jCzyzMcMxAYCNCyZctOo0aNqlbmZcuW0bRp02odG4ea5NWKFWx/6aU0/9e/ePfaa1m6yy45Tre2Ynpu8y1JWSFZeZOUFWqWt0uXLjPNbOeMD5pZ3m7AF8DgDNuXE4pKecc9CqwGGqZtnwKMruz3durUyapr2rRp1T42DjXO+803ZjvsYNa0qdnbb+ciUrmK7rnNoyRlNUtW3iRlNatZXmCGlfO+mu/LVu+T1rYhqR3QmLS2kDTvAYpuax1OKCouVzbaCCZOhObNoXv3sLCUc86lyXfxeAY4SNKGKdv6Ad8D0ys4bnz0tUvZBknNgE7AO7kOWfS22CIMIiwtDb2xli6NO5FzrsDku3jcBfwIPClp/6hdYghws6V035X0kaQRZffNbAbwNDBC0vGSegBjgRXAX/P5BxSN7baDp56CDz+EXr3gxx/jTuScKyB5LR5mthTYD6gPjCOMLL8FGJy2a4Non1THAGOAm4EnCIWja/QzXW3o0iX0vJo+HQYMCEvaOuccMUztYWazgK6V7FOSYdsy4PTo5vLl6KNDu8ell0L79nD99XEncs4VAJ8XylXu4otDARk2LBSQM86IO5FzLmZePFzlJLjtNpg3L6xI2LZtmM7dOVe0fD0PVzUNGsA//hEmU/z97+H11+NO5JyLUaXFQ9LzkraOvj9O0ia1H8sVpCZNYPz4sAJhz57w0UdxJ3LOxaQqZx57A82j7/8G/KLW0rjCt/nmYQZeszAG5Msv407knItBVdo8PgOOlLSMMKJ7y+j7jKLeVK4u+9WvYNw46No1rAXy/PPQuHHcqZxzeVSV4nEtcCdwLmEFwEfK2U/R41VaUtYl3B57wCOPQJ8+oTvv6NFQ3//pnSsWlRYPM7tX0ljgV8CLwBmAn124MPL8L38JPbDOOQduvz30zHLO1XmVFg9JxwETzOwlSVcAT5vZgtqP5hLhrLNgzhy46Sbo0AH++Me4Eznn8qAqDeapjeR/BtrWXhyXSMOGhTXQL7oIqrluinMuWarS5rEUaBN9X9au4dwa9erBgw/CokVhNcLWraFzRUvSO+eSrirFYwrwsKQPovsPSFpe3s5mtmtOkrlkadgwzMK7555w+OHw8suw7bZxp3LO1ZKqFI8TCZMRbg3sBHwKeOd+t66NNw7rgOyxRxgD8uqr0KZN5cc55xKnKr2tSoGbACTtDwwyM1+AyWVWUgITJsA++0CPHvDii7DhhpUe5pxLlqzmtjKzLb1wuErttBM88QS8+y4ccQSsWBF3IudcjmU9MaKkjpKGS3pX0vzo652SOtZGQJdQBx8Md98Nzz0Hp54apjNxztUZWU3JLqkTMA34gbCu+OdAS6AP0F9SFzN7K+cpXTKddFJYB+TKK8MYkMHpC0Y655Iq2/U8bgTeBrpFbSEASGoMTIwer3CVQFdkhgwJBWTIkLCQ1AknxJ3IOZcD2RaPXYG+qYUDQqO6pBuBR3OWzNUNEtxzD8yfD6ecEnpfHXRQ3KmcczWUbZvH90B563lsTLic5dza1lsvNKBvt11oQH/77bgTOedqKNviMQG4TtJeqRuj+9cC43IVzNUxG20EEydCixbQpQu0bUvnrl1D196RI+NO55zLUrbF43zgE2C6pIWS3pG0EJhOGDx4Qa4DujqkTRs44wz45huYPx+ZhUkVBw70AuJcwmTV5mFmS4C9JB0M7AK0BhYCr5vZc7WQz9U1w4evu620FAYNgv7985/HOVct2XbVrW9mq8zsWeDZWsrk6rK5c7Pb7pwrSNletpovaZikbWoljav72rfPbrtzriBlWzzuAo4A/iPpdUkDJW1UC7lcXTV0aOb1zvfdN+9RnHPVl+3cVkPMrCNwAPABcDOwUNLIaNJE5yrWv38Y99GhAyaFM46ddoKHHoIxY+JO55yroqzntgIws+fN7DigFXAW8GtgkqTZkoZI8nm4Xfn694fZs5n+/POht9U//wm77gq//z289lrc6ZxzVVCt4pFiZ2AfwlofS4F/AicDH0k6poY/2xWLxo1h3DjYYgs45BD46KO4EznnKlGdWXU7SBos6WNgKqG77olAGzM7FugA3A3ckNOkrm7bbLOwkBSEGXm/9PXGnCtkWRUPSdOAjwnF4mGgo5kdZGaPmdlPAGa2CniEMNuuc1X3q1+FM5D586FnzzD+wzlXkLI98/gC6A6URI3nc8rZ71/AlpkekLStpKmSSiUtkHSlpPpVDSCpnqQZkkxSzyzzu0K3++7wj3/Am2/C0UfDqlVxJ3LOZZBtb6t+ZvacWcUr+5jZikyFRVILYApgwGHAlYQpTa7IIsbJQNss9ndJc/jhcNtt8PTTcPbZvpCUcwUo2ynZAZDUFtgKaJj+mJlNrODQ04BGQG8z+xaYHI0TGSJpWLStot/bAhgKXALcV53sLiHOPDOMOr/hhrCQ1EUXxZ3IOZci2+lJNgQeAw4s2xR9Tf1oWNElqG7ApLQiMQq4HuhM5bPyXgW8TGiod3XdddfBZ5/BxRdDu3ahK69zriBk2+ZxLdAe2JtQOHoB+wIjCLPq7l7J8VsD76duMLO5QGn0WLkk7UBoqL8wy8wuqerVgwcegH32gQED4IUXYg7knCujSpov1t5Z+gT4E2HFwBXAbmb2ZvTYTUA7M+tbwfErgD+a2a1p2+cBD5nZZRUcO50we+9FkkoIxeoQMxtfzv4DgYEALVu27DRq1Kgq/52pli1bRtOmTat1bBySlLeqWRt89x07nnUWGyxezFu3307plhn7YtS6uvjcFook5U1SVqhZ3i5dusw0s50zPmhmVb4By4G9o++/Aw5KeWw/4OtKjl8BnJth+zzgmgqOOwpYBGwU3S8hXCrrWZXcnTp1suqaNm1atY+NQ5LyZpV19myzVq3M2rUzmz+/1jJVpM4+twUgSXmTlNWsZnmBGVbO+2q2l60+AzaNvv8QSO0quxuVL0O7FGiWYXuL6LF1SFqPMODweqCepOZA2WSMTaJ2GFfXdegQViJcuhS6d4dvK+xb4ZyrZdkWj8lA2QSItwBnSHolGjx4FfBQJce/T1rbhqR2QGPS2kJSNCF0zb2ZUGCWAu9Ej40CfEHsYrHjjmEt9P/8J6yFvmJF3ImcK1pVKh6SGknqQxgk+Iqklmb2MNCH0PawFDiT0IW2Is8AB6WdLfQDvicsZZvJMqBL2q2s281lgC8/V0wOOgjuvRcmTw7L1/oYEOdiUWlXXUkdCQP7SlI2fyupr5k9BTyVxe+7CzgbeFLS9UBHYAhws6V035X0ETDdzE4ys5XAC2mZyrK8a2avZ/H7XV1wwglhDMiQIeFy1pAhcSdyruhU5cxjGLCa0D23MfAbwqWiu7P9ZWa2lNCwXp8wpuMKwuWvwWm7NqDi8SKu2P35z3DiiXDFFTBiRNxpnCs6VRkkuAdwgZm9HN1/T9Kp0dfWZrYwm19oZrOArpXsU1LJ47NZM0DRFSMJ7rorTKJ46qlhOveDD447lXNFoypnHq2BT9K2fUx4826V80TOVdV668Hjj8P228ORR8Jbb8WdyLmiUdXeVt4q6QrThhuGLrwbbww9esDs2XEncq4oVLV4TJL0RdkNKLtUNTV1e/SYc/nVunVYSOqHH8IYkKUZhww553KoKm0e2UyX7lw8tt0WxoyBAw8MU7pPmgQN15n02TmXI5UWDzPz4uGSoXNnePDBMPvu8ceHRaXqZb3SsnOuCqq1nodzBeuoo8I07hddBO3bh/VAnHM558XD1T0XXhgGEd54YyggZ50VdyLn6hwvHq7ukeDWW2HePDjnHGjbFnr1ijuVc3WKXxB2dVP9+jByJOy2Gxx9NLz6atyJnKtTvHi4uqtxYxg7Npx5HHII/O9/cSdyrs7w4uHqts02C2NAJOjWDb7woUjO5YIXD1f3/fKXMH48LFwYzkCWL487kXOJ58XDFYfddoNRo2DGjDAOZOXKuBM5l2hePFzxOPRQuP12GDcOzj7bF5Jyrga8q64rLn/4A8yZA8OGhYWkLr447kTOJZIXD1d8rr02jEK/5BJo1y505XXOZcWLhys+9erB3/4WGtAHDAiz8nbpEncq5xLF2zxccdpgA3jqKdhqqzD6/D//iTuRc4nixcMVr+bNw0JSTZqEdUDmz487kXOJ4cXDFbf27WHChLCAVPfu8O23cSdyLhG8eDj3f/8Ho0fDrFlwxBGwYkXciZwreF48nIOwAuG998LkyXDKKT4GxLlKeG8r58oMGBDWARk8OFzOuvLKuBM5V7C8eDiX6vLLQwG56qpQQE4+Oe5EzhUkLx7OpZJg+PDQ8+q002CLLcJsvM65tXibh3Pp1lsPHnsMdtgBjjwSrr4aSkro3LUrlJSERaacK3JePJzLZMMNQxfehg3Dpaw5c5BZmBdr4EAvIK7oefFwrjytW8P666+7vbQUBg3Kfx7nCogXD+cqsmhR5u1z5+Y3h3MFxouHcxVp3z677c4VibwXD0nbSpoqqVTSAklXSqpfyTG7SPqbpI+i4z6QNFhSw3zldkVq6FBo3HjtbQ0ahO3OFbG8dtWV1AKYAswCDgN+AdxEKGJ/quDQftG+1wMfAjsAV0Vf+9RiZFfs+vcPXwcNwubORU2bwnffwZIl8eZyLmb5HudxGtAI6G1m3wKTJW0EDJE0LNqWyXVmtjjl/guSfgDultTBzObUcm5XzPr3h/79mf7CC+y7996h++6550LbttC7d9zpnItFvi9bdQMmpRWJUYSC0rm8g9IKR5m3o69tchfPuUrUrx+66e6+eygqr7wSdyLnYpHv4rE18H7qBjObC5RGj2VjD2A18HFuojlXRY0awdixYQnbQw+F//0v7kTO5Z0sj7OHSloB/NHMbk3bPg94yMwuq+LPaQX8G5hoZgPK2WcgMBCgZcuWnUaNGlWtzMuWLaNp06bVOjYOScqbpKywbt6G8+ez05lnsqpxY9664w5WtGgRY7q1Jf25LWRJygo1y9ulS5eZZrZzxgfNLG83YAVwbobt84Brqvgz1gdeBD4BWlTlmE6dOll1TZs2rdrHxiFJeZOU1aycvK+/btaokdkuu5gtW5b3TOWpE89tgUpSVrOa5QVmWDnvq/m+bLUUaJZhe4vosQpJEvAQ8Bugu5lVeoxztWrXXeHRR2HmTDjqKFi5Mu5EzuVFvovH+6S1bUhqBzQmrS2kHLcSuvgeZmZV2d+52nfIIfDXv8L48XDWWb6QlCsK+e6q+wzwR0kbmtl30bZ+wPfA9IoOlHQpcCbQ18xeqt2YzmXptNPCpInXXQcdOsAll8SdyLlale8zj7uAH4EnJe0fNWoPAW62lO670UjyESn3jwauIVyymi9p95TbZvn9E5wrx9ChcPTRcOmlPuuuq/PyeuZhZksl7QfcAYwDvgZuIRSQ9FypU5YcGH0dEN1SnQA8kNOgzlVHvXpw//2wcCGccEKYlbdr17hTOVcr8r6SoJnNAir8H2VmJWn3B7Bu0XCu8GywATz5JOy1F/TqBS+9BNtvH3cq53LOZ9V1LteaN4dnnoGmTaF7d5g3L+5EzuWcFw/nakO7djBxInzzDfToAd+WN22bc8nkxcO52vLb38Lo0TBrFvTpAz/9FHci53LGi4dztemAA+C++2DKFDjlFB8D4uqMvDeYO1d0jj8ePvsMLr88rEB41VVxJ3Kuxrx4OJcPgwaFQYRXXx0KyCmnxJ3IuRrx4uFcPkgwfDjMnw+nnw5bbBF6YjmXUN7m4Vy+NGgAjz0WGtL79g2TKTqXUF48nMunpk1hwgTYbLPQhffTT+NO5Fy1ePFwLt9atQqDCH/6Cbp1gyVL4k7kXNa8eDgXh623DkvZzp4Nhx0GP/wQdyLnsuLFw7m47LUXPPwwvPIKHHssrF4ddyLnqsyLh3NxOvJIuPFGeOIJuPDCuNM4V2XeVde5uJ13HsydC7fcEhaSOuecuBM5VykvHs7FTYKbbgqj0M87D9q2DXNhOVfA/LKVc4Wgfn34+99hjz2gf394+eW4EzlXIS8ezhWKRo3g6afDpatDD4UPPog7kXPl8uLhXCHZdNMwBqRBgzAG5PPP407kXEZePJwrNB07wvjxoXD07AnLl8edyLl1ePFwrhDtsgs8+ii89Rb06wcrV8adyLm1ePFwrlD17Al33hnmwjrjDF9IyhUU76rrXCE79dSwDsi114aG9MsuizuRc4AXD+cK39ChYQzIoEFhIaljjok7kXNePJwreBKMGAELFsCJJ0Lr1rDffnGnckXO2zycS4L114cnn4Rf/xp694Z33407kStyXjycS4pmzWDiRNhwwzAGZN68uBO5IubFw7kkadcuFJDvvgtroH/zTdyJXJHy4uFc0uywQ7iE9d57YQLFn36KO5ErQl48nEui/faD+++HqVPh5JN9DIjLO+9t5VxSHXtsWAfkT38KXXivvjruRK6I5P3MQ9K2kqZKKpW0QNKVkupX4bhmkv4maamkbySNlLRJPjI7V7AuuwxOOSWMBbnnnrjTuCKS1zMPSS2AKcAs4DDgF8BNhCL2p0oOfwzYCjgZWA1cD4wB9q6luM4VPilMYTJ/Ppx+OmyxBfToEXcqVwTyfeZxGtAI6G1mk83sLuAK4HxJG5V3kKQ9gAOB481stJk9BRwD7CVp/3wEd65gNWgQJlHccUfo1Qtat6Zz165QUgIjR8adrmIjR0JJSTLyJikr1HrefBePbsAkM/s2ZdsoQkHpXMlxn5vZi2UbzOwN4NPoMeeKW9OmYfT5ypWwaBEyC3NiDRxYuG9yI0eGfHPmFH7eJGWFvOTNd4P51sDzqRvMbK6k0uixcRUc936G7e9Fjznnhg1bt9dVaSkcdxycf364xAXZf63OMVX5+tFH6041X1oKAwbANddU+ufm1f/+l5ysUH7eQYPCMsc5kO/i0QL4OsP2pdFj1TmuY6YDJA0EBgK0bNmSF154IYuYayxbtqzax8YhSXmTlBUKP2/nuXNRhu22ejULdt8dIHwKhTVFJq3YVPR4TY7N9Phm77+fOe/KlXy52WYZHonPZrNmJSYrVJB37lym5+o1bGZ5uwErgHMzbJ8HXFPBcZOBMRm2/x14pbLf26lTJ6uuadOmVfvYOCQpb5KymiUgb4cOZuGtee1bhw5xJ8ssSXmTlNUsZ3mBGVbO+2q+2zyWAs0ybG8RPZbr45wrHkOHQuPGa29r3DhsL0RJypukrJCXvPkuHu+T1kYhqR3QmMxtGuUeFymvLcS54tO/fxjr0aEDJoXFo+65J2fXuHMuSXmTlBXykjffxeMZ4CBJG6Zs6wd8D0yv5LhWkvYq2yBpZ0J7xzO1EdS5ROrfH2bPZvrzz8Ps2YX75lYmSXmTlBVqPW++i8ddwI/Ak5L2jxq1hwA3W0r3XUkfSRpRdt/MXgWeAx6S1FvS4cBI4CUzm5LPP8A551yei4eZLQX2A+oTuuVeAdwCDE7btUG0T6p+hLOT+4GHgJlAr9rM65xzLrO8T4xoZrOArpXsU5Jh29fACdHNOedcjHxKduecc1nz4uGccy5rsiJYREbSl8Ccah6+KbA4h3FqW5LyJikrJCtvkrJCsvImKSvULG8HM8s4hL4oikdNSJphZjvHnaOqkpQ3SVkhWXmTlBWSlTdJWaH28vplK+ecc1nz4uGccy5rXjwql7S1PZOUN0lZIVl5k5QVkpU3SVmhlvJ6m4dzzrms+ZmHc865rHnxcM45lzUvHs4557LmxcM551zWvHg455zLWt5n1XW5Ea3A2B0Q8LiZLZHUFrgQ+AUwG7jHzN6NLyVIuhiYGHeOqpLUCGhgZt+lbNsMOBPYFlgN/Au408y+iSWkcwXAu+pGJImwPkgPYBtgY8IbxSLgNeABM/tffAnXkLQrMBloAqwEvgIOAiYCq4D/AtsBrYD9zeyfMUVF0mrACMsFPwKMMrOP48pTGUkTgQ/N7Jzo/h6E1SpXE9aQEdAJ+Anoamb/jTHrjkAjM3slZdvBwKWsKXTvAENS9ykU0f+5Q4CdCK+RGYQPGgX9piRpI8JcUV3N7KW488DPmboC6wMTzGx59KHnDMKKq58QPkwuyNnvLPB/p7yInuSJhDeFRYQ3hi0IL+hnCE/+r4GrzOyquHKWkTSZcNbYC1hOWFDrcMKb2xFmtkLSBsAYoKGZdYkpalnxuB7YHjiAkPstQiF5zMzmx5UtE0mLgZPM7Ono/muE5/jwsrMRSc2AscAPZnZQjFlfA8aZ2dDo/onAfcA04HlCodsP2BvoU/Y3xZT1FcLz+l50vwVhddBOwLJot6aED2oHpZ75xUHSHyp4uBFwA/AX4EMAM7szH7kykfRLYCrQLtr0KXAg4QNmc+BjwvvX90AnM5uXk19sZkV/A/5BeBFsn7KtDfAsMDq635nwIj+xAPIuAbql3N+c8CnzwLT9egCLY866Gtg1+n5jYGD0Ql8Z3V6Itm0c9/MaZSwF9km5/1P685ry3C6POeu3qdmAj4DbM+x3F/BOobwOovsjCGfMB6dsOxhYCtxSAK+D1YSz+NXl3FIfWxVz1scIZ5i/jP6PPRy9n70CbBjts2m0z925+r3eYB50Ay6xlOvyFk7vTgMOl9TazKYD1wDnxJQxnWX4Pv00sqBOK83sKzO7x8z2A9oCFxBOs+8CFkqaEGvA4D9A6pna54T/kOk2IRSaOK1Ou98BeCLDfk8QPnkWkkOBK83s2bIN0fdDgd6xpVrjaeAL4CSgvpnVK7sRXg8C9o22pS+ZnW97AUPN7CMz+wr4E6Hd80aLzuDMbDFwK2u/tmvEi0dQj/BJIt0qwoukWXT/dWCrfIWqwAzgQklNJdUDLgPmA6dLqg8gqQHwB8KbYcExs0Vm9hcz+x2wJWEd+zYxxwK4DrhE0onRczgUuEHSAZLWl7RB1K5wLeETX5z+CfRPuf9fINPU27sQXh+FpDnhdZxuJqGtLlZm1gs4Hvgj8KakPVMfjidVuVoQLreXKfu3Tl/D6BPCh7ac8N5WwWTgakn/NrNP4OdrsrcR/lHKGsqbAoXQw2YQIfNSwqWfUkJj2ePAh5LKGszbEC4FFDQzm0N4076uALI8Kekswqe0W4APCB8eniW8aSjadSzhjSVOlwEvRx8gbic0lD8oaWPC5UARXhfnApfElDFVH0llxW0pkGmRoU0Jl+NiZ2bPSdqB8PxNkPQsoTdjrO0xGXxBOOssswq4m3DWnGpzcpjdG8yBqIvrJMJZxRzCde4tgR+B35vZM9F+wwgra/WLK2uZKHNPwgeA0Wa2UFIr4CLCJYo5wH1m9laMMZE0GLjXctjLIx8kbQL0A3YlfBIW4Q3vPWC8mc2MMd7PJP0fMBzYjbWLW9n3SwmXh/4SS8BI1HEi3QNmdmLafncD25rZ3vlJVjXR/61hhEtqdxMKShczezHOXACSxgBfpT+XGfa7HdjazA7Iye/14hFEl3v6Ar8FGhIaHx+JriE6V9AkbUMoIOmF7hUzWxFntmxIOgX42MyejztLJlHX7VsIH9B6WAF0gZbUEmhsZp9Wst/5hI4TU3Pye7141D2S6plZpk96BUNSQ0Kj3mrgo0J8g4vaPDqSMubHzObGm8q5wuAN5mkk/UZSH0knR7c+kn4Td650knpLGiNpoqRDom39JM0GVkiaE32Ki5WkY6LxB2X3G0i6jtBN89+EBv2vJBXCNXkAJHWSNJZwffg94GXgVeBTSfMlXSmpcawh6xBF4s6RiaRG6f/Wkv4vel/oFFeughBn/+RCugEnEtoJMvXtXkWY7uOEuHNGWftGuV4idCksBU4htNWMIIwq/UeU+6CYs84CTk+5f1OU93JgT0I3wyGEAUyXFcBzeyChrWsGoWv2EMJA0Z+izBcQejX9C2hRAHl7EsbNvBe9FvbJsM9uxD8W4UCiMQcp2w4nDBhdCayInvMecT+nUbZmwFNRrpXAvUB94MG094WXgE3jzlvFv6lPLl8Hsf9BhXADzopeJH8ljMbdNHqh1I++3wu4I3pTOaMA8r4J3JVyv3+U7aa0/f4GTIk5aynQOeX+F8A5Gfa7EJhTAM/tTODBcl4jswln6w2jN707Y856QPQG9nL0+pwZ3b+J6JJ0tF8hFI9VrD1IsFf0BvxK9G9/AeHsbiUZBmXGkPc2whQkZwHHRR8YRgOfRYVwM8L4sPnA8LjzVvFvymnx8DYPQNInhDfjYZXsdxFwmpl1zE+ycnN8C/Q2synR/WaEBtL9LaWhMbqcdbeZxTZ+QtJC4EwzGx3d/5FwNvRC2n4HAGPNrFH+U66V43vgUDObnLa9BWFk/2/M7D1JxwHXm1nrOHJGmV4izMN1Qsq2EwlvfJMJPQV/kLQboeE8tsFsUW+r3c3sjej+W8B8Mzskbb+JQBMz6xxDzNQcswkD7+6N7u9IKM4nmNmDKfudQjhj3jKWoCHD/VXctQNhYGNOXgfe5hG0Bt6own5vUAADmAjdMFNfAGVzA32dtt8ywmCsOI0lDGhcP7o/Bfh9hv1+T/h0F7cvCD3u0v2W8LyXjfOZw5rBo3HZDvh76gYzu58wlc7uwPPRmI9CtB2hy2u6ewgTJcZtM9aM74JoDivCPFGpPiLzeJV8Op5wNrR9JbcO5f2A6vBBgsE7wCmSXrRyeilFDXqnEBp54zaHMGvqJAAzWxV1IXwvbb9fAAvznC3dpYSR0P+RdB8wDrhe0nasGcjWBdiRMMNq3O4BrpLUlFDofiKM0B4ETLM141U6AnH3vPqBMLPyWsxsZjQiehLhstCQPOcqT+pljm/IPGBtOYXxofZTQhGeHt3fm3CZ7XeEdo4yexL/6+BD4A0zO66inSQdATyaq1/qxSO4gDCCeJakJwnTh38dPdYM2JpwjbYthTFi+0nSphkws9cz7HcUa7/Q887MvpK0O+HN93zCKFeAPaLbT4RLLHub2ZvxpFzDzIZGl1guAf5ctpnQAeHclF1XEBrU4/RvwnX3sekPmNknUQGZCDyQ51zlmSRpZfR9M8IHhulp+2xN/B94IMy39hdJ2xMKXV/CB6HLow8W7xDOkM4D4p5p+zVCUatM6iDSGvM2j4ikXxBGZx/MmqmNy3xG6HFzgxXwWhTpJLUHvjazgpjuAUBSCWsPZPvYCnOMx3qEM7eGwCeF9ByWkXQqYYqSHa2cwaySmhB6De1vYVK/WEQzDaT70MweSdvvhWh7IXQzP5twOXU9wmwNd0n6PaFNqWxizHuAi+N8DUddhvc0s9sq2W9TQptdesGu3u/14rGuqF938+ju12YW9+ypzrkCEV3C3tTMvow7S5y8eNQx0Sn1W0D/QrgMpAQu66qELPHrXJy8eKSI3jTaEC6lLM7w+KZAdzN7KO/h1s7RvYKHmxAaxS4hmo7dzCbmI1cmStCyrpCsJX6rKpr36kgzuzLmHHlfKjWXojOO1GVzZxL+jtjfRBVmK+5D+P/0gJm9L+m3wBWs+cBzh5lNytkvjXvgSiHcgA0I05mvim4rCCO1m6XtF/tgqyhHklY5WwwclnL/NcKI6A1TtjUjNJxOKoDndjJhGdfmhGvddwDzCKO310t5vTxD6H0V++u3Cn9TTgeHVTPDLwm9BMtelx8T3tQ+IRToNwlTsX8OtC2A5+wVYJuU+y2ijKujnN+yZpDjhnHljLIdRPjwtSh6Xr8l9GBcGuX7a/T/bhVhOeXc/N64/5EK4UboVfM1oSvuzsDZ0Yv4Q+BXKfsVSvGYQeiRcgKh73bqbYfoRd23bFvMWROzrGuUI0lL/Lav4u20uF+3xLRUag3yJmbZXMIMA48TVjyE0IliKTAibb+Hgddy9nvj/kcqhBuha+6ZadtaAS8CXwJ7RNsKpXiIsO73F4QpE7ZMeaxZ9MJfZ46jmLK+AQxOuf8ZcFSG/Y4DviyAvEvS3iA2i57PA9L2614AxaPsLLOyWyGcgS4A+qbc7xDl6p223wnA/wrgdZBePL4Ezs2wX+zT6hC6Eu+fcr9FlL9r2n4HEjoA5eT3+jiPoB1pg//MbJGk/QjVeoqk/hRG/3MsvBLukfQYcDXwbrTQy9XxJsvoOmCkpM+Ah1izrOsSwqWqskGChbCsK6xZ4vclwllT6hK/z1sYkFkoS/x+BzwP3FfJfnsRuqHHKZalUnOoOYW7bO73rD1YtOz79Kl+GhMGluaEF49gAfArwpnGzyz03T5K0l8Ip4WxNpSnM7OvgTMl3UPoe/4hcD0FtMayJWtZV0jWEr9vENrlJlS0U7R2StxiWSq1hpKybO7LwJ8lfRhluZEwm/XF0awZ30Xz311EKHY54b2t+HlisY5mtm8F+1xK+NRsFuMEcxWRdBRhqcy2hAnQYl8is4wSsqwrJGqJ38uBgWaWPqg1fb99gCvMrEt+kmXMMIYYlkqtLiVo2VxJvyRMpVP2OphNOJt/gjBifw5QQvgw1MXM/pWT3+vF4+dubv2Aa62CZWclHU249n1CefvELbqk0gRYZmar4s7jHMS3VGptU4EsmxuN79qT0ENwqpl9Hw12Ppk1H3geMbN5OfudXjycc85lqxBmr3S1RNK9kkbEnaMqkpQVkpfXuVzzBvMsSLoXqGdmJ8WdpYq6kJwPCEnKCgnKK2kK4SrDfnFnqUySskKy8uY6qxeP7CTmDQPAzH4Zd4aqSlJWSFxekZzXbZKyQrLy5jSrt3nUYVEXzc3NLO7FaiqVpKyQvLzO5VpSKmZBkNQwWiMjKXoQVkRLgiRlhQTllbReUl63ScoKycqb66xePLKTmDcMVxwknSHpY0nfS3pH0rEZdtuJAnjdJikrJCtvHFm9zSOBJFW1T3mmEbF5laSskKy80aDQ2wlL5L5NWIr0AUmHAceYWc6moqipJGWFZOWNK6sXD5L1hhHZhzDNx6xK9iuEaSmSlBWSlfdC4EYz+3neqmg+tpHANEk9zWxJbOnWlqSskKy8sWT1BnNA0kqq9oaxBbBb3NOTSHoHeN/M+lWy3xHAo3HmTVLWKEdi8kr6DjjEzF5I215CWG+kPmH+rc2AVzxr1SUpb1xZvc0j+C/wHzM7sqIbcHPcQSOvAbtXYb/UiQfjkqSskKy83xAm5luLmc0mXLpYDLwK7JLfWBklKSskK28sWf3Mg58nNzvYzDpUsl8fwprWsRZdSb8AfmNmYyvZrxGhO2n6tNd5k6SsUY7E5JX0NPCdmR1TzuONCJPjdSPmCT2TlDXKk5i8cWX14kGy3jCcKyPpSOA8oGd5E3pKqg8MJ0zouWU+86XlSEzWKEti8saV1YuHc865rHmbh3POuax58XDOOZc1Lx6u6EgaIGmmpO8kLZX0tqRa6UknaStJQyQ1r8K+QyRZym2BpNFRm1xlxz4gKdMa287VCi8erqgoLCd8HzAJ6A0cBzwNHFpLv3IrYDDQvIr7fwPsEd0uBP4PmCqpSSXHXQUMqFZC56rBR5i7YnMmcLeZXZaybZykK+IKlGalmb0Wff+apLnAP4HuwOPpO0tqZGbfm9nH+QzpnJ95uGLTHFiUvtFSuh1KKokuGx0t6eHo8tYXkganHyepq6TXJf0g6XNJdyqsJ42kfYFx0a6fRj9zdpZ5Z0ZfS6KfOVvSTZIulzQP+Dbavs5lK0kdJP1D0mJJpZL+LenolMcbShom6TNJP0YT6nXPMp8rUn7m4YrNW8BZ0Sf68ZXM+XMDMB44gjDn1WBJi83srwCSfgM8C0wG+gDtgOuAjoTpIN4imneIcIlsIfBjlnlLoq+pBe9owqwIf6Cc/8OSNieMKi6NMnwGbBdlLPMEsCvhstrHQF9grKSdzexfWeZ0RcaLhys2ZwBjgAcAk/QeMJowsdy3afv+18xOjb6fFL0hXyZpuJmtBi4H5gCHmtkqAElfAY9K2sPMXpX0QXT829F0EZWSVPb/siNwJ/AdMCVtt56VzJZ6HtAM6GRmC6NtU1N+x36EJQb2NbPp0ebnJG0FDAKOrEpWV7z8spUrKmb2b2AbQgP5nYT5qS4HZpRdbkrxVNr9J4E2QNvo/q7AU2WFIzIaWAnsVc2ImwArotsHhALSL6UAAEytwjTbXYFn045LtT/hbOZlSQ3KboQCs3M1s7si4mceruiY2Y+EtohxAJJOIvTAOgn4S8quX6QdWna/NTA3+vp52s9eJWkJsHE1431DeGM3wpv7Alt3GojP1zlqXZsAb1bw+KZAK0KRSrcqwzbn1uLFwxU9MxshaRiwddpDm5dzf2HK17X2ieYQ2gTIOMdQFaw0s8rGa1RlTqElhOJWnq+A+cDhVczl3Fr8spUrKlG7Rfq2zQjtA+mf6Hul3S9r9J4X3X8d6BUVjNR9GgAvRfd/ir7me/GoqcBBklpW8HgrYJmZzUi/5S+mSyo/83DF5t1oCuvnCJehOhB6I5UCD6bt+5touv7RhN5WJwHnRI3lAFcTlv0cI2k4oS3kemCSmb0a7VPWYH6qpFFAqZm9Wzt/2lpuIQyA/KekoYTeVtsATcxsGKGH2CRgsqTrCb23NiIMSmxoZpfmIaNLMC8erthcCRwG3EZol1gEvEJolP40bd+LgJ6E4vEDYRT3HWUPmtl/JXUDriE0pn9LWEf6opR95ki6EDgbOItw1lJSG39YKjP7UtKewDDgVmAD4EPg2uhxk9QbuAw4F2hPuJT1L8J62M5VyKdkdy5NtHznp4SlPcfHHMe5guRtHs4557LmxcM551zW/LKVc865rPmZh3POuax58XDOOZc1Lx7OOeey5sXDOedc1rx4OOecy9r/Az15xUncHan3AAAAAElFTkSuQmCC\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, strike_price - x)\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:37:33.823468Z",
     "start_time": "2020-07-13T23:37:33.818301Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "exact expected value:\t0.1709\n",
      "exact delta value:   \t-0.8193\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": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:37:35.232511Z",
     "start_time": "2020-07-13T23:37:35.227091Z"
    }
   },
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "problem = EstimationProblem(\n",
    "    state_preparation=european_put,\n",
    "    objective_qubits=[num_uncertainty_qubits],\n",
    "    post_processing=european_put_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": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:38:02.357252Z",
     "start_time": "2020-07-13T23:37:38.158364Z"
    }
   },
   "outputs": [],
   "source": [
    "result = ae.estimate(problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:38:02.364454Z",
     "start_time": "2020-07-13T23:38:02.359842Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact value:        \t0.1709\n",
      "Estimated value:    \t0.1770\n",
      "Confidence interval:\t[0.1720, 0.1820]\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 a comparator circuit and an ancilla qubit to identify the cases where $S_T \\leq K$.\n",
    "However, since we are only interested in the (negative) 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": 10,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:38:02.785334Z",
     "start_time": "2020-07-13T23:38:02.780118Z"
    }
   },
   "outputs": [],
   "source": [
    "# setup piecewise linear objective fcuntion\n",
    "breakpoints = [low, strike_price]\n",
    "slopes = [0, 0]\n",
    "offsets = [1, 0]\n",
    "f_min = 0\n",
    "f_max = 1\n",
    "\n",
    "european_put_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",
    ")\n",
    "\n",
    "# construct circuit for payoff function\n",
    "european_put_delta = european_put_delta_objective.compose(uncertainty_model, front=True)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:38:02.794225Z",
     "start_time": "2020-07-13T23:38:02.787236Z"
    }
   },
   "outputs": [],
   "source": [
    "# set target precision and confidence level\n",
    "epsilon = 0.01\n",
    "alpha = 0.05\n",
    "\n",
    "problem = EstimationProblem(\n",
    "    state_preparation=european_put_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": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:38:25.893471Z",
     "start_time": "2020-07-13T23:38:02.795889Z"
    }
   },
   "outputs": [],
   "source": [
    "result_delta = ae_delta.estimate(problem)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {
    "ExecuteTime": {
     "end_time": "2020-07-13T23:38:25.900397Z",
     "start_time": "2020-07-13T23:38:25.895473Z"
    }
   },
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Exact delta:    \t-0.8193\n",
      "Estimated value: \t-0.8197\n",
      "Confidence interval: \t[-0.8236, -0.8158]\n"
     ]
    }
   ],
   "source": [
    "conf_int = -np.array(result_delta.confidence_interval)[::-1]\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": "2020-07-13T23:38:26.340110Z",
     "start_time": "2020-07-13T23:38:26.189871Z"
    }
   },
   "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_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_finance</code></td><td>0.4.0</td></tr><tr><td><code>qiskit_ibm_provider</code></td><td>0.6.1</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:14:14 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
}