{
 "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",
    "- <a href=\"https://arxiv.org/abs/1806.06893\">Quantum Risk Analysis. Woerner, Egger. 2018.</a>\n",
    "- <a href=\"https://arxiv.org/abs/1905.02666\">Option Pricing using Quantum Computers. Stamatopoulos et al. 2019.</a>"
   ]
  },
  {
   "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",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "\n",
    "from qiskit import Aer\n",
    "from qiskit.aqua.algorithms import IterativeAmplitudeEstimation\n",
    "from qiskit.circuit.library import LogNormalDistribution, LinearAmplitudeFunction"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Uncertainty Model\n",
    "\n",
    "We construct a circuit factory 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 factory 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(num_uncertainty_qubits, mu=mu, sigma=sigma**2, bounds=(low, high))"
   ]
  },
  {
   "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": "iVBORw0KGgoAAAANSUhEUgAAAZwAAAEyCAYAAADOV2anAAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3debgcVZnH8e9PEAgEMLKERSSACiJRMEGMoCSArDpBEKKiYxCNjAKOIgwiQgAdAWVxYBQjakRGoyKiIIhhSRDZ40JYomwhBhQEAjGEQELe+eOcC51K973dt/tWNcnv8zz93FtVp6re7tu33646myICMzOzgfaKqgMwM7OVgxOOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpKk84kraVdI2khZIekXSKpFX62OdNkn6Tyz8naY6kCyRtXCg3WVLUeWwzsM/KzMyKVq3y5JKGAFcDdwNjga2AM0mJ8IRedl0XeBC4EHgE2AI4CRghaceIWFJTdhZwaGH/2c3Et/7668ewYcOaKTognnnmGdZaa63Kzt9It8YFjq0/ujUucGz9UXVcM2bMeDwiNqi7MSIqewBfAOYB69SsOxZYWLuuyWO9GwjgrTXrJgO39ze+ESNGRJWuu+66Ss/fSLfGFeHY+qNb44pwbP1RdVy9feZWfUttH+CqiJhfs24KMAjYtcVjPZF/rtaJwMzMrLOqTjjbkG55vSgi5pCucPqsZ5H0CkmrSdoaOA24Dbi1UGxbSfNzXc8NklpNZGZm1gGKCsdSk7QYOCYizimsnwtcGBHH97H/b4C98uIMYN+IeKxm+2eA50l1RBsARwMjgF0iopiYevaZAEwAGDp06IgpU6b056l1xIIFCxg8eHBl52+kW+MCx9Yf3RoXOLb+qDquMWPGzIiIkXU3NrrXVsYDWAz8Z531c4H/bmL/1wM7AR8mXSnNANbopfyapMYGlzYTn+tw6uvWuCIcW390a1wRjq0/qo6LLq7DmUdqcVY0JG/rVUTcGxG3RMRFpCudHYAP9VJ+IXAF8Nb+hWtmZv1VdcKZRaGuRtJmpCuRWXX3aCAiHgKeBLbsq2h+mJlZiapOOFcCe0lau2bdOOBZYHorB8oNB9Yj3TJrVGYQsB/p1puZmZWo0o6fwPnAUcAlkk4nXZ1MBM6KmqbSku4DpkfEYXn568AS4BbgKeCNpP4795OaVSNpXeBy4CLgPmB94LPAJsBBJTw3MzOrUWnCiYh5knYHzgMuIyWPs0lJp9aqQO1wN7cDR5Jak60BzAF+Dnw1Ip7JZZ4D/kkasWBDYBFwE7BrRNw+EM/HzMwaq/oKh4i4G9itjzLDCstTyFcyveyzCDig3fjMAIYd9+u2j3H08CWMb/M4s0/br+04zKpSdR2OmZmtJJxwzMysFE44ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVJxxJ20q6RtJCSY9IOkXSKn3s8yZJv8nln5M0R9IFkjauU3aspJmSFkm6W9K4gXs2ZmbWyKpVnlzSEOBq4G5gLLAVcCYpEZ7Qy67rAg8CFwKPAFsAJwEjJO0YEUvy8XcBfg58EzgK2Bf4saR5EfHbAXlSZmZWV6UJBzgcGAQcEBHzgamS1gEmSjojr1tORNwI3FizapqkucBvgTcDf8jrvwRcHxFH5eXrJL0JODGXNTOzklR9S20f4KpCYplCSkK7tnisJ/LP1QAkrQ6MAX5aKDcFGCVp3dbDNTOz/qo64WwDzKpdERFzgIV5W68kvULSapK2Bk4DbgNuzZu3Al5ZPD5wD+l5v6G90M3MrBWKiOpOLi0GjomIcwrr5wIXRsTxfez/G2CvvDgD2DciHsvbdgZuAHaIiD/V7PM64F5gr3r1OJImABMAhg4dOmLKlCn9fXptW7BgAYMHD67s/I10a1wwcLHNfPjpto8xdBA8+mx7xxi+aecvzFfGv2cndGtsVcc1ZsyYGRExst62qutw2nUk8Grg9aRGBldK2jkiFvX3gBExCZgEMHLkyBg9enQn4uyXadOmUeX5G+nWuGDgYht/3K/bPsbRw5dw5sz2/uVmHzK67TiKVsa/Zyd0a2zdGhdUn3DmkVqcFQ3J23oVEffmX2+R9DtSy7UPAd+r2b94/CE15zYzs5JUXYczi0JdjaTNgDVZvu6lVxHxEPAksGVedT+wuHj8vLwU+Gs/4jUzs36qOuFcCewlae2adeOAZ4HprRwoNxxYj3SVQ0Q8B1wHHFQoOg64KSLavylvZmZNq/qW2vmkDpmXSDqddHUyETirtqm0pPuA6RFxWF7+OrAEuAV4CngjcCzpqqa2lv9UUh+dc4BLSR0/9wX2HtinZWZmRZUmnIiYJ2l34DzgMlLyOJuUdGqtCtQOd3M7qcHABGANYA5pRIGvRsQzNce/QdL7gS8D/0Gu4/EoA7aiGNahxgztNoqYfdp+bcdhK76qr3CIiLuB3fooM6ywPIVlr2R62/dS0tWNmZlVqOo6HDMzW0k44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZlcIJx8zMSuGEY2ZmpXDCMTOzUjjhmJlZKZxwzMysFE44ZmZWCiccMzMrRcujRUsaDrwN2Ig0NcCTpNkzb4wIT9tsZmZ1NZVwJG1Jmk/mEGAoaYrmp4DngFeRpoReKmk6cAHwk4hYOiARm5nZy1Kft9QkXQDcBWwPnALsAKwRERtExGsiYjCwIfBeYCZwBnCPpF0GLmwzM3u5aeYK51lgm4h4qFGBiHgcuBK4UtLngIOATTsTopmZrQj6vMKJiCN7SzZ1yi+NiJ9ExE+aKS9pW0nXSFoo6RFJp0hapY99dpT0fUn35f3+IukkSWsUyk2UFHUeezf7fMzMrDPammJa0nbAroCA6RExs8X9hwBXA3cDY4GtgDNJifCEXnYdl8ueDtwLvBk4Nf88sFD2aaCYYO5pJU4zM2tfvxOOpP8AvgJcA6wFfE3S0RHxzRYOczgwCDggIuYDUyWtA0yUdEZeV89p+TZej2mSFgHflrR54YpsSUTc3EJMZmY2AJppNLBmg03/BYyKiIMiYl/gCOCLLZ5/H+CqQmKZQkpCuzbaqZBsevwx/9ykxRjMzKwEzXT8/KukQ+qsF6l5dI/+NIPeBphVuyIi5gAL87ZWjMox3F9Y/ypJj0taLOmPkg7oR5xmZtYmRUTvBaR3AecAzwNHRcStef2nSc2kryH1w9kdODYizm365NJi4JiIOKewfi5wYUQc3+RxNgLuAK6IiPE16z9MarL9R2Bt4JPAvsCBEXFJg2NNACYADB06dMSUKVOafTodt2DBAgYPHlzZ+Rvp1rhg4GKb+fDTbR9j6CB49Nn2jjF803WXWe7WuDplZXyvtavquMaMGTMjIkbW29ZnwgGQJOAwUsX8VOC/IuLvkt7CS7e+ro+IP7USWCcSjqTVSA0PXgOM6G20g/w8bgQGRcT2fR175MiRcfvtt/dVbMBMmzaN0aNHV3b+Rro1Lhi42IYd9+u2j3H08CWcObOtdjrMPm2/ZZa7Na5OWRnfa+2qOi5JDRNOU2OpRXIBsDXwKHCnpC8CsyLif/KjpWSTzQPqfTUakrf1KieQC4E3Afv2NbROpOx6CfDmvppem5lZZ7U0eGdEzI+IY4CdSOOpzZL0/jbOP4tCXY2kzUi36GbV3WNZ55CaU4+NiGbKA0R+mJlZiZpqpSbpy5JuyZXuk4BFETGWVNdxkqTp+fZaq64E9pK0ds26caTRDab3EdcXSC3jPhwRNzRzsnxFdCDw54h4oR/xmplZPzVz4/a7wLakPjcLSUlmqqRtI2JqTjSfyusujYgJLZz/fOAo4BJJpwNbAhOBs2qbSku6j9Sx9LC8/CHgv4HJwMOS3l5zzPsj4p+53HTg56SrpbWAT5CuzvZvIUYzM+uAZhLOPsBBETEVQNLvgSdIPf3vy6NCnyfpR8BJrZw8IuZJ2h04D7iMNAL12aSkU4yzts5lz/xzfH7UOpSUiADuA/4T2JjUZPoPwH4RcWUrcZqZWfuaSTizgI9ImgEsIjUtfgaYW1soIp4EPtNqABFxN7BbH2WGFZbHs3yiqbffYa3GY2ZmA6OZhPNR0hXD46TK9tmkK55FAxeWmZmtaPpMOBHxF2CUpLWA1Tyrp5mZ9UfTvb0i4hnSrTQzM7OWNdMs+iOtdpKU9DpJ7+x/WGZmtqJppuPn54D7JZ3aW18bSetJOkTSZcCfSC3DzMzMgObqcHaQNA44EviipAWkCcweB54DXgVsAbyWNBzNRcDhEfHwgEVtZmYvO03V4eTpon8iaStgD+CtwEakzpSPAtcDvwemRcTiAYrVzMxexloaIjYi7mf5+WbMzMz61NLgnWZmZv3lhGNmZqVwwjEzs1I44ZiZWSlaSjiS3ivJScrMzFrWavK4FJgr6XRJbxyIgMzMbMXUasLZCvgOcDBwp6SbJH1C0jqdD83MzFYkLSWciJgdESdFxBbAu0kTnJ0N/F3SDyWNGYggzczs5a/f9TERcW1EfAR4AzADOAS4WtIDkj4rqaVOpWZmtmLrd8KRtKukycBfgO2A/yVN/XwxcDJwYScCNDOzFUOrrdQ2l3SipPuBa4HNgAnAxhFxZERcExHHkmYJHdvkMbeVdI2khZIekXRKX9MhSNpR0vcl3Zf3+4ukkyStUafszpJukbRI0oOSjmrlOZuZWWe0etvrAeAR0pTT34uIBxuUuwu4ta+DSRoCXA3cTUpQWwFnkhLhCb3sOi6XPR24F3gzcGr+eWDN8V8HXAVcDnwBeBtwlqSFEXFBX/GZmVnntJpw3gNcFRFLeysUEX8FmmlAcDgwCDggIuYDU3OLt4mSzsjr6jktIh6vWZ4maRHwbUmbR8RDef0xpAT54YhYAlwr6bXASZK+GxHRRIxmZtYBrdbh7EialmA5kjaWdGKLx9uHlMBqE8sUUhLatdFOhWTT44/55yaF41+Sk03t8V9DqncyM7OStJpwTiJ9WNezSd7eim2AWbUrImIOsDBva8UoYCl5+gRJa5HqmGYVyt1Tc24zMyuJWrmrJGkpsFNE3FZn21jguxGxfgvHWwwcExHnFNbPBS6MiOObPM5GwB3AFRExPq/bFJgLvC8iLq0puyqwGPhkREyqc6wJpIYQDB06dMSUKVOafTodt2DBAgYPHlzZ+Rvp1rhg4GKb+fDTbR9j6CB49Nn2jjF803WXWe7WuDplZXyvtavquMaMGTMjIkbW29ZnHY6kj5JanQEE8C1JxbqVNYDhwG/bCbQ/JK0G/BRYAHy23ePlJDQJYOTIkTF69Oh2D9lv06ZNo8rzN9KtccHAxTb+uF+3fYyjhy/hzJntdU+bfcjoZZa7Na5OWRnfa+3q1riguUYDC4En8u8CngaeLJR5HrgS+GaL558H1PtqNCRv65Ukkfr7vAnYOSJq93kq/ywef0jNuc3MrCR9JpyI+BnwMwBJ3wdO6aU5dKtmUahLkbQZsCbL173Ucw6pOfW7I6JYF/SMpL8Vj1+z3MzxzcysQ1odS+3QDiYbSFdFe0lau2bdOOBZYHpvO0r6AnAEqcnzDb0c/32FjqTjgL8Bd/Y7ajMza1nVc9ucDzwHXCJpj1xhPxE4q7apdB5R4Ls1yx8C/pt0O+1hSW+veWxQc/yvkVrV/VDSGEnHAp8kXaW5D46ZWYmaaTRwKzA+Iu6WdBup4UBDEfG2Zk8eEfMk7Q6cB1xGqnc5m5R0inHWXqXsmX+Oz49ah5JGQiAi7pO0N3AW6WrnH8DRHmXAzKx8zTQauIt0i6vn945eGUTE3cBufZQZVlgez/KJptG+N5CGtDEzswo102jg0Jrfxw9oNGZmtsKqug7HzMxWEs3U4fRZb1OrlTocMzNbeTRbh+MWXWZm1pZm6nDGlxCHmZmt4FyHY2Zmpai0H46Zma08Ku+HY2ZmKwf3wzEzs1K0PAlGnn9mPKn3/sbA34FbgB9ExPMdjc7MzFYYLTUakPRG4F7gf4HtgBfyz/8F7pO0bccjNDOzFUKrVziTSBOwvTMi5vSslPRa4HLS6M/v6lx4Zma2omg14YwEPlibbAAiYo6kk4AfdSwyW+kM69B0ye1Ouzz7tP3ajsPMltdqP5zZwBoNtq0BzGmwzczMVnKtJpzjgC9L2ql2paS3A6cC/9WpwMzMbMXSn8E71wFulPQY8BiwYX48ARwPXDoAcZqZ2ctcfwbvvGuAYjEzsxVY5YN35qbU5wKjSFNMXwCcHBEv9LLPasBXgLeTGjKsERGqU24y8NE6h3hjRMxqP3ozM2tWyx0/O0nSEOBq4G5gLLAVcCapbumEXnZdE/g4cCtwI71PUT0LOLSwbnb/IjYzs/6qNOEAhwODgAMiYj4wVdI6wERJZ+R1y4mIpyS9OiJC0hH0nnCeiYibOx+6mZm1ouXpCSSNk3S1pDmSHis+WjzcPsBVhcQyhZSEdu1tx4jwIKJmZi8jrQ5t8yHgB8B9wGuAX5FGGHgFMB84r8Xzb0O65fWi3Kl0Yd7WCdtKmi/pOUk3SOo1kZmZ2cBQKxcKkv4IXAycBiwGRkbEHyStDUwFLo6Ir7dwvMXAMRFxTmH9XODCiDi+iWMcAZzboNHAZ4DnSXVEGwBHAyOAXSLi1gbHmwBMABg6dOiIKVOmNPt0Om7BggUMHjy4svM3MlBxzXz46baPMXQQPPps3+V6M3zTdZdb162xdWtcndKt/wPQvbFVHdeYMWNmRMTIettarcN5PfD7iHhB0gukPjlExL8knQ6cDTSdcAZaRHyjdlnSFaRm3ccD+zfYZxJpzDhGjhwZo0ePHuAoG5s2bRpVnr+RgYqr3SFpIA1tc+bM9qomZx8yerl13Rpbt8bVKd36PwDdG1u3xgWt1+HMB1bPvz8MvLFmm4D1WjzePKDeV6MheVtHRcRC4ArgrZ0+tpmZ9a7VrzW3AW8GriLV35woaQnpttWJQKutwWZRqKuRtBmp2fNA9ZMJPGupmVnpWk04XwU2z7+fmH//FulK6Tbgky0e70rgGElrR8S/8rpxpCmtp7d4rD5JGgTsB8zo9LHNzKx3LSWc3J/l5vz7U8BYSasDqzfqM9OH84GjgEtyHdCWwETgrNrjSboPmB4Rh9Ws2wdYC9g+L78/b7otIh6StC6pBd1FpFZ16wOfBTYBDupHrGZm1oaOTTEtqeUppiNinqTdSc2pLyMNbXM2KekU41ylsO5bvHS1BfCz/PNQYDLwHPBP0ogFGwKLgJuAXSPi9lbiNDOz9rWUcPIU078hXSXMII0WvR3w78CXJO0dEXe3csxcvreRAoiIYc2sK2xfBBzQSixmZjZwPMW0mZmVotVm0SOBE+tNMQ2cBOzYqcDMzGzF4immzcysFK3eUjsOOFPSgxFxS8/KmimmP9/J4Mzs5WtYh0ZBaHc0hdmn7dd2HNYZnmLazMxK4SmmzcysFJVPMW1mZiuHfg0RK2kTYBTwatKttJsj4pFOBmZmZiuWVjt+rgKcC3yCZXv+vyBpEnBkRCztYHxmZraCaLVZ9MnAx0iNA4aRpoIelpc/xvJD0piZmQGt31L7d+CEwqyec4CvSQrSQJwndio4MzNbcbR6hbMhcEeDbXfk7WZmZstpNeH8FfhAg20fAP7SXjhmZraiavWW2peBKXmwzouBR0lXNQcBY2icjMzMbCXX6gRsP5X0FKnxwDeAVwKLSVMV7B0RUzsfopmZrQiaTjiSXkmadO3OiBgl6RWkWTQfd1NoMzPrSyt1OC8A1wLbAETE0oh4zMnGzMya0XTCyYnlXmCjgQvHzMxWVK22UvsicKKk4Z0KQNK2kq6RtFDSI5JOySMa9LbPapK+Jul3kp7NfYAalR0raaakRZLuljSuU7GbmVnzWm2ldgKwHvAnSQ+TWqkt82EfEW9r9mCShgBXA3cDY4GtgDNJifCEXnZdE/g4cCtwI7Bbg+PvAvwc+CapU+q+wI8lzYuI3zYbp5mZta/VhHMXcGcHz384aXicAyJiPjBV0jrAREln5HXLiYinJL06IkLSETRIOMCXgOsj4qi8fJ2kN5FGQ3DCMTMrUavNosd3+Pz7AFcVEssU4HRgV+CyXmJpeBsNQNLqpL5BRxU2TQG+L2ndiHi6X1GbmVnLmqrDkTRI0oGSjpb0IUlDO3T+bYBZtSsiYg6wMG9rx1akfkKzCuvvIT3vN7R5fDMza4H6uFBA0pakepZhNavnAwe3Ww8iaTFwTEScU1g/F7gwIo5v4hhHAOdGhArrdwZuAHaIiD/VrH8dqbXdXvXilzQBmAAwdOjQEVOmTGn9iXXIggULGDx4cGXnb2Sg4pr5cPsXnEMHwaPPtneM4Zuuu9y6bo2tW+OC7o6tE1a2/89mjRkzZkZEjKy3rZlbamcAS4F3kkYU2IJUCf/t/PsKJSImAZMARo4cGaNHj64slmnTplHl+RsZqLjGH/frto9x9PAlnDmzX/MKvmj2IaOXW9etsXVrXNDdsXXCyvb/2QnN3FIbRZqS4PcRsSgi7gE+CbxW0sZtnn8eUO/rx5C8rd1jU+f4QwrbzcysBM0knI2BBwrr7gdE+51AZ1Goq5G0GanZc7HupVX3k8Z5K9YFbUO6Yvtrm8c3M7MWNNvxs/eKnv67EthL0to168YBzwLT2zlwRDwHXEcaybrWOOAmt1AzMytXszdHr5K0pM76a4rrI6KVSdjOJzVbvkTS6cCWpGmqz6ptKi3pPmB6RBxWs24fYC1g+7z8/rzptoh4KP9+KjBN0jnApaSOn/sCe7cQo5mZdUAzCefkgTp5RMyTtDtwHqnPzVPA2aSkU2tVoDjczbeAzWuWf5Z/HgpMzse/ISeiLwP/ATwIfMijDJiZla/PhBMRA5Zw8vHvpvFIAT1lhjWzrsG+l5KubszMrEKtDt5pZmbWL044ZmZWCiccMzMrhROOmZmVwgnHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVJxxJ20q6RtJCSY9IOkVScTrpevutK+n7kuZJelrS/0lar1BmsqSo89hm4J6RmZnV0+cU0wNJ0hDgauBuYCywFXAmKRGe0MfuPwXeAHwcWAqcTppK+p2FcrOAQwvrZrcTt5mZta7ShAMcDgwCDoiI+cBUSesAEyWdkdctR9IoYE9g14i4Pq97GLhF0h4RcXVN8Wci4uaBfRpmZtaXqm+p7QNcVUgsU0hJaNc+9nu0J9kARMStwIN5m5mZdZmqE842pFteL4qIOcDCvK3p/bJ76uy3raT5kp6TdIOk3hKZmZkNEEVEdSeXFgPHRMQ5hfVzgQsj4vgG+00l3Srbv7D+ImDLiHhHXv4M8DypjmgD4GhgBLBLviKqd+wJwASAoUOHjpgyZUobz7A9CxYsYPDgwZWdv5GBimvmw0+3fYyhg+DRZ9s7xvBN111uXbfG1q1xQXfH1gkr2/9ns8aMGTMjIkbW21Z1Hc6Aiohv1C5LugK4Czge2L/BPpOASQAjR46M0aNHD3CUjU2bNo0qz9/IQMU1/rhft32Mo4cv4cyZ7b2tZx8yerl13Rpbt8YF3R1bJ6xs/5+dUPUttXlAva8fQ/K2ju4XEQuBK4C3thCjmZl1QNUJZxaFOhdJmwFrUr+OpuF+WaO6nVqRH2ZmVqKqE86VwF6S1q5ZNw54Fpjex34bSdqlZ4WkkcCWeVtdkgYB+wEz2gnazMxaV3XCOR94DrhE0h65wn4icFZtU2lJ90n6bs9yRNwE/Ba4UNIBkvYH/g+4oacPTh6J4HeSPilpd0njgOuATYD/LusJmplZUmmjgYiYJ2l34DzgMuAp4GxS0qm1KlAc7mZcLvs9UuK8HDiqZvtzwD9JIxZsCCwCbiJ1Fr29o0/EzMz6VHkrtYi4G9itjzLD6qx7ijRkTXHYmp7ti4ADOhCima1ghnWoBV07LfFmn7Zf2zG83FR9S83MzFYSTjhmZlYKJxwzMyuFE46ZmZXCCcfMzErhhGNmZqVwwjEzs1I44ZiZWSmccMzMrBSVjzRg5eqGHtawcvayNlvZ+QrHzMxK4YRjZmalcMIxM7NSOOGYmVkpnHDMzKwUTjhmZlYKJxwzMyuFE46ZmZWi8o6fkrYFzgVGAU8BFwAnR8QLfey3LnAOsD8pcV4OHBURTxTKjQW+DLweeCAf+yedfh5mZu1a0TtmV3qFI2kIcDUQwFjgFOBo4OQmdv8pMBr4ODAe2BG4tHD8XYCfA9cB+wC/Bn4sac+OPAEzM2ta1Vc4hwODgAMiYj4wVdI6wERJZ+R1y5E0CtgT2DUirs/rHgZukbRHRFydi34JuD4ijsrL10l6E3Ai8NuBe1pmZlZUdR3OPsBVhcQyhZSEdu1jv0d7kg1ARNwKPJi3IWl1YAzpSqjWFGBUviVnZmYlqTrhbAPMql0REXOAhXlb0/tl99TstxXwyjrl7iE97zf0I14zM+snRUR1J5cWA8dExDmF9XOBCyPi+Ab7TQWeiYj9C+svAraMiHdI2hm4AdghIv5UU+Z1wL3AXhGx3G01SROACXlxa+Av/X6C7VsfeLzC8zfSrXGBY+uPbo0LHFt/VB3X5hGxQb0NVdfhdJ2ImARMqjoOAEm3R8TIquMo6ta4wLH1R7fGBY6tP7o1Lqj+lto8oF5dypC8rZ39en4Wyw0pbDczsxJUnXBmUairkbQZsCb162ga7pfV1u3cDyyuU24bYCnw137Ea2Zm/VR1wrkS2EvS2jXrxgHPAtP72G+j3M8GAEkjgS3zNiLiOVL/m4MK+44DboqIp9sPf8B1xa29Oro1LnBs/dGtcYFj649ujavyRgNDgLuBO4HTSQnjLOCciDihptx9wPSIOKxm3VWk0QM+T7piOR14LCLeWVNmF2AacB6pU+i+ufze9RoMmJnZwKn0Cici5gG7A6sAl5FGGDgbOKlQdNVcptY40lXQ94ALgRnA+wrHvwF4P7AHcBXwb8CHnGzMzMpX6RWOmZmtPKquwzEzs5WEE1I1XdsAABkPSURBVI6ZmZXCCcfMzErhkQa6gCSRGjzsB7wReDWp5d0/gJuByRFRSb+h3C9qX0DAzyLiCUmvIbX22wqYDUyKiJklxvRfwBVlnrNZkgYBq0bEv2rWbQAcAWxL+rv+Cfjmy6RpfiXy/8R7gbeSpi+5nfQ374pK5zyq/ePAbrlxUlUx7AasBvw6Ip7J77VPk1r8PkD633ykivjqcaOBiuU3yBXACFKCeR7YlPRPdiXpjbM1cGpEnFpybG8DpgJrAUuAJ4G9crwvAHcB2wEbAXtExO9Kimsp6fWZBfwImBIR95dx7r5IugK4NyI+k5dHkf6OS0ktKUX6Wz9P+rC6q6S4dgAGRcSNNev2Br7AS4nwz8DE2jIlxXYjcFhE3JOXh5CmDxkBLMjFBpO+fO1Vm8wHOK5P9bJ5EPA14BuksRmJiG+WERe8OCbkNcBmedWDpClbpgKvInV835rUp3FERMwtK7ZeRYQfFT6AH5PesMNr1m0C/Ab4eV7elfSP97GSY5tK6jz7KtLI2+cBc4FfAq/MZVYnfaBeV2JcS4GvkmZ5fY6U/G4DPgtsWvHf83FgbM3yzaQPhrVr1q1LatJ/VYlx3Qx8sWb5Y/l1vAb4InBC/lsvqY2/xL/n22qWv0v6crN3zbq9ScNRnV1yXC/kn/UetdteKPk1+ynpC8LrSHdEfpg/R27sea+RBvH8M/DtMmPrNe6qA1jZH6RptQ+ss35YfkNvnJePB/5ccmxPAPvULG+Y/7n2LJTbD3i8xLhe/IDK/2wT8gfnkvyYlte9uoK/50LgXTXLzxdfr5rX7JkS45pfGwdwH3BunXLnV/A+KyacfwL/Wafc54GHSozrF8DfgUPJd4Nqtr0qx/2usuIpnP8R4OCa5c1zPAcUyh0K/LWKGOs93Gigeq8gJZaiF0i3X3oGH72FaubwiTq/F+/DVnZfNiKejIhJEbE78BrSFOWrkT44/y6p/UniW3MnaeK/Ho+SkmLReqTkVJalheXNgYvrlLuYdCumSq8i1dkUzSDdvi1FRLwP+ChwDHBbnvLkxc1lxdHAENIt+B4P558PFco9QPq/6ApOONWbCnxZ0pY9K/I97P8hvaF6GgsMBsquZL4d+LykwZJeQbrKehj4D0mr5FhXBT5F+qCtVET8IyK+ERHvALYgjVixSclhnAYcJ+lj+bX5CvA1Se+WtJqk1XPdyVdZfjbagfQ74JCa5buAekPY78hLH15lOlDSp3K9yTyg3nwq65Ou1EoTaVSSN5NmCv61pCm50UzVHiN9aejxAvBt0hecWhsCpdR5NaXqS6yV/UH69nEXaWTr+0hjyz1LutVWezvrDOAnJcc2kvTPvzjH9ATwFlISfIA0HNGDpHqUMSXGtcwtmG57AB8nfTA+Ddyaf3+BdLvvhfz4BbBmiTENz3H8EHgbaSr2x0gJ8d2kCufTgEXUuZ1Vwt+z+PhenXLfBn5X4d91I9IwWguAM/PfsapbapfWe43qlDsXmFrVa1Z8uJVaF8hXCweTPszXICWeH0XEk5UGBuRvc+8hNaH/eUT8XdJGwLGkWy8PARdExB9KjOkk4DvRRc09iyStRxrv722kDyqRkvc9wOURMaOCmLYHvgXsRLolpLyp5/d5wCkR8Y2yY2uGpE8A90fEtRXHMYo05uPWwH5Rcqu+HMNQ0heWB/so9zlSndw15UTWOyccs5WMpDeSkk4xEd4YEYurjM1WbE44XUTSm0gTxNXOSjorSuqr0SpJr4iIYmV0ZSStQeqMuhS4r+oPz1yHsyU1HXkjYk6VMb3c5A6gRIUfVLkzryJiYc267ckdn6u4Wn25cqOBLpArmB8C7gB+RppAaVL+/Q5JsyUdWlFsB0i6VNIVkt6b142TNBtYLOmhfKujzJg+LOljNcurSjqN1HfjDlIDhiclHVdmXDXxjJD0K1Jl7T3A74GbgAclPSzpFElrVhFbN5K0Z2ESRiTtL+kPpPrD5yXdLmm/kuNaV9IvSHVf8yV9R9Iqkn4A/IH0/3mrpBskrV9mbM2SdKCkeq1gK+GEUzFJR5IqQy8HRpNalbwyPzYkdfq8HDhf0qdLju1gUjPZ9Un/+D/JyeWHpH4vR5E6mp0vaa8SQzue1OG0x+k5lq8C7yK9ZmcCJ0k6vsS4kLQn6TXZhHSf/1RSr/kXgImkCQYPBG7MrRHLjO09kq6RdI+kX0p6V50yO1XwAXUlaUinnhjeB1xCasBwXH4sBn6ZX9+ynAq8E/gcqaPsO0gtC3cjdUQdSqrf3CKXtT74llrFJD0AnB8RZ/RR7ljg8IjYsrdynSTpNmBGRByelw8hTXh3XkQcXVPu+8BmEbFHSXEtJLXgm56XHwO+UqzslvR54MiI2LzOYQYqthnAnRHx0cL6I0l9hLYk9RO6Ebg5InobPqWTcb2bNHrFzcAfgVHA9sA5wOd7bllJ2olUl1Oc8HAgY1sKvD0ibs3LfwAejoj3FspdAawVEbuWFNds0vvqO3l5B1JfoEMj4gc15T4BHB8RW5QRVz7n95osujkwusy/Z298hVO9jUlNZ/tyKyV2esu2ZtnOgZeTrryKnSkvIY3HVZanSVddPdYlDeFR9GfSVWKZtgUuqrP+IuC1wNYRsYj0Qf++OuUGyknAhRGxc0QcEREjgE8AnwQuyfVf3WI70lV/0STSYJ5l2YCX+sFBHjONNE5Zrfuo329oIH2U1JR9eB+P0r5sNcMJp3p/Bj6RO1bWlStOP0GqnyhTsOzU3j0DKT5VKLeA1Du8LL8idUhdLS9fDXywTrkPkvo4lekxUvP2oreQXs+ezrsP8dIoEmXYjkIijIjvkW4/vh24VlK9ERHKUnur5Wnqd1Z8hnI/sx4kvT493klq/PGOQrmdgbIbg9wLXBsRO/b2IN2O7BqenqB6R5Nuddwt6RLSCMg9H+jrklqtvY/UQXTvkmN7iPSN/SqAiHgh90G4p1BuK9KYU2X5Aqnn/J2SLiB1QD1d0nakcdREGl5mB9IQ92WaBJwqaTApET5P6r3/RdIApz19h7ak3A+pRaRRv5cRETPykC1XkW7zTSwxplpXSVqSf1+X9LebXiizDeW+z84HviFpOCkJHkx6730p/33/TLri+izl1+HczPKJr57a/laVc8KpWET8PjexPJY09MhmhSJ/I1Wqfi3KH4L/EgrjMEXELXXKfQAobU6QiHhS0ttJH+Kf46XbZqPy43nSkEHvjIjbyoorx/aVXCdxHHBiz2rSqOD/WVN0MfDfJYZ2B2l0gV8VN0TEAznpXAFMLjGmHifXWfdYnXUHkka0LkVEnJfvPHyQ1DDg2Ig4X9Jc0tBTPePhnQ98vay4snNJLeX6Mp1lx/arlBsNdJncXLbn9tRTtW3/u5Wk15JiLXWcq5rzD2PZToz3d0EfnFeSrvzWAB6o6rWpieeTpNZ9OzQawULSWqQhd/aICN9u70W+zb1+RPyz6lheTpxwzMysFL6l1iWUpnLehPTt/PE629cH9o2IC0sPro58D/sPwCFl37ZSl0/jrC6clrvb6WUyXXK+sqmd+noGKd7Sv7lLGkm6zSjSNPSzJL2FdIuy5312XkRcVXZsjfgKp2KSVie1Hjogr1pKGpH2c7UflhX1j9i3l81rAT8h1VXcCRARV5QUV1dO45xj6cppuZulNM7aQRFxSonn7MrpktWlU1/nWPYiNZZ5ktR6bwNgLKne9R5SX6sRpAYrB0bEpWXF1qsyhqT2o9fhw08ktUr7BGk6gKNIc1rcC7y+ptxOlD+NbVdOsUuXTuOcz9uV03K3EP+BFbzPunK6ZLp06ut83t+ThtZZJS8fn+P4bqHcD0kdjCt/b0V4iunKH6Rm0EcU1m0EXE+aandUXldFwrmdl6bY3bzweHP+hzy4Z12JcXXlNM75nN06Lfdrm3wcXsH7rCunS66TcLpi6ut8zqdJV8g9y0NyvLsVyu1JatBTWmy9PVyHU73NKHTojIh/SNqd9O3k6jykTJn9D3rsSLryOp3U7+XzkeffkNTTafEfEVGc1nag9UzjfH1e7pZpnHt047Tcs5s8p5os10kvl+mSu2Lq6+xZlu1X1fP7oEK5NUl9sLqCE071HgFez0sfngBEatb7AUnfIF06l95YINJXpEmSfgp8GZgp6dz8e5VOA/5P0t9Ir0vPNM5PkG6j9XT8LHsaZ3hpWu4bSMmudlruayN1nq1iWu5/AdcCF/RRbhdSn7AydfN0yQfmynnooqmvSbfUTpR0bz7310mzBf+XpOsj4l/5S+GxpITYHaq+xFrZH6TBMKf1UeYLlFxP0iCON5N68j8CfIZqp9jtummcc1zdOi33VNJQKH2Vq6IOpyunS6aLp74m1XfNrnmv30+6JdrzvzCTlJznAduXGVtvD1/hVO+bwDhJr44GHfIi4qtK8+W8u9zQlovjDmC0pA8AZ1DhkBkRcUGeq6RnGucn6YJpnCPi9jwUSnFa7nfx0rTcV1LytNykK+gJTZT7J4Wr7RJ8knTrpy8PkpJTKaL5zq+3k1psliYi7stDOe1MapxyTUQ8K2k06cvY1qRb8j+Kklr1NcPNoq1f8m2htYAFEdE1EzyZWfdywjEzs1J4vKSXiTy97XerjqOebo2tW+OC7o6tW0m6WtI1VcdR1K1xQffF5jqcl48xdO8XhG6NrVvjgi6NTdLVpDsfu1cdSx2iC18zujcu6LLYfEvNzF6Uvw2/IiK6Zkh7W3F0Teaz3klaI08D0HW6NbZujQu6N7aI2L1bk42kV3bja9atcUH3xeaE8/KxH6lZaDfq1ti6NS7o0tiq+oCS9GlJ90t6VtKfJX2kTrG3UvJr1q1xdXtsjTjhmK0kuvUDKvfrOpc0COuXSJ0YJ0u6WNIaZcbycoir22PrjRsNVEzStU0WrTekxoDq1ti6NS7o3thqPqB+TBq6/h2kD6ixwIcjosrxtj4PfD0iXhxSJ48l+H/AdZLeExFPOK6XTWwNudFAxSQtAf5CGgepN5sCO0W58+F0ZWzdGhd0b2ySbicNbVPvA+pB4D2RJoqrYt6lfwHvjYhphfXDSKMyrEKaBmCDMmPr1ri6Pbbe+AqnencBsyJiXG+FJL2fkofPoHtj69a4oHtj25r0rfhFEXGNpLeTPqBukrR3ifHUepo0AOYyImK2pHcAvwZuAk51XC/q5tgach1O9W4G3t5EuaD8scu6NbZujQu6N7aGH1Ck22uPkz6gdiwxph4zgP3rbYiIecDupPHK/qfMoOjeuKC7Y2vICad6ZwBHNlHuCmCLAY6lqFtj69a4oHtj6+YPqIuALSXVm9OIiHgW+DfS1ApzHBfQ3bE15Docs5WApIOAz5LqauqOSi5pFeBbwLsjouxEbSsBJxwzMyuFb6mZmVkpnHDMzKwUTjhmZlYKJxwzMyuFE471StJ4STMk/UvSPEl/lHTWAJ3rYEnjmyg3UVLUPB6R9HNJWzV5nsm5533lmn3OuWzP8763wfZ78/aJAxVDi8dd5nXu9HkkvULSEfk9+ayk+ZLukvQ/kvrVx0nJnyR9tMH2ybk3f71t58mT6vXKCccakvQFUjv+q4ADgH8Hfklq3z8QDgbGN1n2aWBUfnwe2B64RtJaTex7agvnGWitPGeARcAWkkbWrpS0IzAsbx/oGJpVfJ07fZ6fAF8GLiG9Jz9K6t/0juh/89uDgVcDP+rHvl8HDpH0un6ee4XnoW2sN0cA346I42vWXSbp5KoCqrEkIm7Ov98saQ7wO2Bf4GfFwrmPySoR8XxE3F9inJ32DPAH4AOkjpo9PgBcC4yoIqgeZb3OkvYB3g/sGxFX1mz6RX+vbrKjgB9GxOKac61KSp4fATYBPijpfuDkiHhxeKI8rMwNwH8AR7cRwwrLVzjWm1cB/yiurP322HPbRNL+kmZJWiTpBknbFvfLt1RmSnpO0t8kfSX/MyNpMnAgsGvNrbKJLcQ6I/8cVieuu0jf/Heq3VaI7V2SrpO0QNLTkqZJ2qFm+zslTZe0UNITkr4jae3eApI0StKvJP1d0jP5Vs0hta9dP5/zFODgng/W/PPgvL5jMeTX4OLC8UbnMtvVvpZ9vc6NziNpX0lLJW1ROM8Wef3YBq/BrvnncqNz9/fqJl+ZvAO4uLDpM8CxpFEYrgA+BnwPWK/OYX5OusrxZ2sdvsKx3vwBODJfPVzey3DnmwNnkebleBY4GbhK0ut7hr2XtCfpFsiFwDHAm0nfGtcDDs+/v5aU5D6Vjzu3hViH5Z//KKw7Azglr687z4uk0cBU4DrSbZlngJ1JIzr/UdLOwNXApaRv1esBpwFD8nIjmwO/B84nfRDvDHxf0tKI+DH9f86XkEYE2IV0VfdO0qjAlwBfKymGWsPo+3VudJ6/A4+QXveJNeXHA4+RBqGs55n882uSzoyIh1qMuZ7d83H/XFi/K2mk7TPyF6nf5zHo6rkRGAoMr3Mciwg//Kj7ICWFB0gDTS4ljYR8CrBOTZnJefs7atZtDiwBDq9ZdzNwXeH4xwIvAK/JyxcD05qIayJpsMlV8+MNpGQxH9i4ENf2dfafDNxes3wT6faUGpzvd3Vi3y0ff7smX0vlWL9N+vDqWd/Uc6593vn3XwL/m3//JnBp/v1xYGInYgCmARcX1o2ufd4tvs6NzvNlUpJSTZyzSfO9NHotNgLuyOcO4E7geGBwG+/3ScBtddZ/G/hbPudkYFgvx1g1v/c/0d84VuSHL/usoYi4A3gjqUL2m6QPgi8Bt0saXFP0sYi4sWa/h0i3uN4GL97XfyvL1638hHRbd1Q/wlsPWJwffwG2BMZFxN9ryjwcEX/q7SC5kcFOwA8if2IUtq+Z4/uppFV7HsAN+dwN60wkDVFqMfVQTawTSAmyXVOA90tanXSVtdzttBJi6NHn69yH75G+pIzOy2Py8vcb7RAR/wB2APYiXe29CvgKcKOk1QAkfTTfQvxTvo07K/8+Q9Ir6xx2I1LCLvoK6crnQdL/wufzVW+9uJYAT+VjWYETjvUqIp6LiMsi4oiI2Bb4OPB64LCaYo/V2fUxYOP8+/rAK4FHC2V6luuOeNuHp0lD6Y8EXkP61nlloUzxfPUMISXSv/eyfRVSwl1c83iO9Jw26+XYk4FxpNtce+Z4vwd0YgrgXwGDSR+GawGXVRBDj2Ze54Yi4gHS1dShedWhwK0RcVcf+70QEb+NiE+Rbtd9n3Qra1Te/oOI2J70ZWcJsHNEbB8RI6KmUUCNNUh/1+J55uTjvo90xb8LcIMadw94js6+visM1+FYSyLiu5LOALapWb1hnaIbkm7BQfrWuLhOuaH5Z93Ri/uwJCL66kvTTOXxPNLtwo0bbH8qH2ciqcK46JF6OynNK/8e4NMRcX7N+o58yYuIZyRdThoB+mcR8UyxTAdiWASsVlg3pF44TR6vNxcA31Fqin8ALbbyioilkn5LSlbFD/vXA/Oi7ymXn6TBlUlOUL9Rmqp7Immqh7MlnZMTUq1X0b/39ArPVzjWkKTlEomkDYB1WfZb7YZKswz2lHkt6VvlrZC+iZJusR1UONzBpA/7m/Ly85T8zTB/UN8C/HtPq686228Gto6I2+s86iYcYHXS/9eL35hzq7ZiH6Z2nvO3SFc25zfY3m4Mc1n2iwWkq6T+6u25XpK3TyHFXPcWIYCkoQ02/RuwkPT3rPUWmqvA/wt15iiq974Abss/X10ouwGwJvDXJs630vEVjvVmpqRfAr8l3SLbnNTJciHwg5pyjwMXSTqBl1qpPUa6ndPjJFLLte+TPkyGk1oufScielpFzQLGStqf9GH3SC8f6J10HKkV2pWSJpHu148iVXhfTmrccI2kpaSK73+RbuHsB3wxIpb7cImIpyXdBpwoaT4psR5HuhW4Tk3Rfj/nSPPZT+tle7sx/AI4TNLZpNZiY4B2pqFu+FwjYpGk/wM+Dfw4Ip7q5Tg/lfQv4KekxgUbAocAY0mV9cV930JqYNCX35Neqw0i4p81638k6Y/A9aTblyNIV5YPA/cUjjGSdMV3I7a8qlst+NG9D9I//29Jt40Wkf65fwRsU1NmMqmF1wGkb3XPkf5xl2u9RapLmEn6JjuXVP+was329Ukfck+Sb2M1iGsiubVWL7FPpqaFVF/bSE1frycl06dIrd62r9m+E/AbUku4Z4C7SU3B1+0lhtcB1+Tyc0iJa5nYm33OLTzvZVqptRsD8AVSC61/kWaZ/DeWb6XW1Ovc13MF9sjr9+jjOX4s/y3m5vfSk6SEOLpB+cuADzTxfl8NeAL4SGH9+/L5/kFK2vNJiX6HOsf4BoUWjX689PAEbNaW3KFvu4gY2VdZs97kusGDgS0jYmkHjzsH2Csiilcj9cp+A3hdROzXYPtkUqKcXWfbKsBDwHERcVFbQa+gfEvNzColaWtgW9KQMCd3ONkMIXWKbbZO5WvAXyW9IercKu3DQaRbyg3rn1Z2bjRgZlX7NulW7RWk4WM6JiLmRcSgSA1Xmik/l3TLrlGrxUtJt1zrEXBYpL44VodvqZmZWSl8hWNmZqVwwjEzs1I44ZiZWSmccMzMrBROOGZmVgonHDMzK4UTjpmZleL/AW46+1xJ9j1lAAAAAElFTkSuQmCC\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",
    "<a href=\"https://arxiv.org/abs/1806.06893\">Quantum Risk Analysis. Woerner, Egger. 2018.</a>"
   ]
  },
  {
   "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+CAYAAAB1DJw3AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4yLjIsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy+WH4yJAAAgAElEQVR4nO3dd5hU9dnG8e8NqDRFbIhSVpIYNWpexRoLgpViARSMWLChxl5iIwYsWLBGjdgwlmCwYKEpAiLGLmiMCWpsgDQVxAJroTzvH78zMgyzZXZn58zZeT7XNdfunDln995hmGfO+TWZGc4551wuGsQdwDnnXPJ48XDOOZczLx7OOedy5sXDOedczrx4OOecy5kXD+eccznz4uGKmqTBkiztNk/SKEm/iClPM0kjJS2K8vSPtp8s6VNJyyW9UMGx+2T8Lanb8kL+DVGWLaPndv2M7f2jTM0LncklS6O4AzhXDd8AB0XfdwCuBCZL+o2ZLS1wltOAg4FjgbnAx5I2BYYBtwOPAYur+Bn9gE/S7scx2GpLYBBwP/B12vZxwO5AeQyZXIJ48XBJsNzMXou+f03SbOCfQDfCm3UhbQV8YGajUhsk7Qk0BO4zs39X42f828z+U1cBa8PMvgS+jDuHK35+2col0fToa1l0Gel2SR9IKo8uHf1V0nqpnSU9mu1SUnTZ5nNJa0X3N5L0QHRJqlzSC5J2Stt/JnAisEPaJafBhEIG8E76paxcRb/v8YxtqUtd20b3y6L7fSTdJekbSXMkXS6pQcax20saI+lrSUskvSFpf0n7AGOi3T6Nft7M6Jg1LltV9byknhtJN0g6N8qzOLq8t9plMVd/ePFwSVQWfV0ANCV86h8IdAUuA7qw+hnJcGBvSVukNkgScBzwdzNbFm1+CjgQuADoS/j/MUXSL6PHewLjgfcJl3Z2B+4FTo8e7xdtG1dF/oaSGqXdavL/cCiwBDgc+Dvw5+j71N+3FfAy0Bo4Ncr+JNAWeCv6GwF6RZl7VvK7qnpeUvoA+wIDgIuAHsDVNfjbXBKYmd/8VrQ3YDCwkHCJtRHhWv0U4FugdZb9GwF7ENoR2kXbGgCzgcvT9usS7bNtdP+g6H6ntH2aES7h3JW27X5gWsbv3Cf9Z1Xyt6T2y7xdFT3+AvB4ZT+bUDgNeDBjv38BI9Pu/wOYAzSpIEuP6OeUZWzvH21vnuPzMhP4GGiUtu0WYEHcryG/1c3N2zxcEmwILEu7Pxvoa2bzASQdA5wH/IrwxpayJTDbzFZK+htwrKTBFt7Z+hOKQKrtYRfgCzObmjrYzJZKGgvsmee/50jCG23KvBr8jOcy7s8A2qXd70I4q/q+Bj87XS7PyxQzS+85NgPYRNJatursztUTXjxcEnwD7Ef4BLwAmBcVACT1BB4k9Ha6FPiKcKnmSaBx2s/4G+GSVmdJbwK9WXXphuiYL7L87s+BDfL5xwD/tdo3mH+dcf8nVv97NwTm1/J3QG7PS7ZMAtZh9eLv6gEvHi4JlpvZtAoeOwJ43cz+kNogqVPmTmY2U9IkwhnHFoRLWf9I22U+sEmWn9+KUJAK4Qdg7YxtLWv4sxYR3vhrqxieF1eEvMHcJV0T4MeMbf0q2Hc44YzjD8BTZpb+Sfl1wiWWvVMbJDUFugMv5S9upeYQugKnO6CGP2sy0EdS4woe/yn6WtHjKcXwvLgi5MXDJd1EQk+qgZL2k3QTocdPNk8RPt3vSLiM9TMzmwC8Ajwi6ThJPQg9q5oA19dZ+tU9CfxK0s3R3zKEVYMjc3U50AJ4UVLf6Of9UdIJ0eMfRF9PkbSrpO2y/ZAieV5cEfLi4ZLuLuBG4GzgCaA9cFS2Hc3sR+AZ4DNgUpZdDiMUo1sIXX0FdDGzj/IfO2u+cYR2m8MJhaQ94e+qyc/6gNCgvZDQnfjJ6OfOih6fRWjz6UXo0jsm+08CYn5eXHFS1O7oXL0nqRHhzfM+M7ss7jzOJZk3mLt6T9LawG8JZyQbEs5WnHO14MXDlYLNgDcIXU5PMbM5MedxLvH8spVzzrmceYO5c865nJXEZauNNtrIysrKanTs0qVLadasWdU7Fokk5U1SVkhW3iRlhWTlTVJWqF3e6dOnLzSzjbM+GPfkWoW4dezY0WpqypQpNT42DknKm6SsZsnKm6SsZsnKm6SsZrXLS8YkoOk3v2zlnHMuZ148nHPO5cyLh3POuZx58XDOOZczLx7OOedy5sWjIiNGQFkZnbp0gbKycN855xxQIuM8cjZiBAwYAOXlCGDWrHAfoF9FS0U451zp8DOPbAYOhPLy1beVl4ftzjnnvHhkNXt2btudc67EePHIpl273LY751yJ8eKRzZAh0LTpmtsPP7zwWZxzrgh58cimXz+4+25o3x6ToG3b0ONq2DCYNi3udM45FzsvHhXp1w9mzmTq88+Hto5XX4VNNoHu3eGTT+JO55xzsfLiUV2bbgrPPAPLlkHXrrBoUdyJnHMuNl48crHVVjB6dBj3ccgh8P33cSdyzrlYePHI1Z57wt//Hi5jHXMMrFgRdyLnnCs4Lx41cfjhcOONMGoUXHBB3Gmcc67gfHqSmjr33HD56pZbwviPc8+NO5FzzhWMF4/auPFGmDMHzj8/dOf1cSDOuRLhxaM2GjaEhx6C+fPh6KNDj6w994w7lXPO1Tlv86itJk1CD6z27eHQQ+H99+NO5Jxzdc6LRz5suGEYA9KoURgDsmBB3Imcc65OefHIlw4dYNw4+OIL6NEDliyJO5FzztUZLx75tNNO8Mgj8Pbb0LcvLF8edyLnnKsTXjzyrUcPuOMOGD8e/vAHMIs7kXPO5Z33tqoLp5wSJlO8+urQkO4rEDrn6hkvHnXlqqtCAfnTn8IYkGOPjTuRc87ljRePuiLB8OFhDMiJJ8Jmm8F++8Wdyjnn8qLgbR6StpE0WVK5pHmSrpDUsBrH7STpOUlfRbdJknYtROYaW3vtMP/V1ltDr17wzjtxJ3LOubwoaPGQ1BKYBBhwKHAFcD5weRXHtY2OawQcE90aARMlta/LzLXWokVoPF9vPejWDT77LO5EzjlXa4U+8zgVaAL0MrOJZnYnoXCcJ2m9So7rDqwL9DSzcWY2DugJNAO61XXoWmvTJgwiXLIkFJCvv447kXPO1Uqhi0dXYIKZfZu2bSShoHSq5Li1gOXA0rRtS6JtynfIOrHddvDEE/DBB+ES1k8/xZ3IOedqrNDFYytgtcmfzGw2UB49VpFR0T43StpE0ibAzcBi4LE6ypp/++4bGtGnTIETTvAxIM65xCp0b6uWQLZrNoujx7Iys3mSOgNjgbOizfOBA83sy7ynrEvHHBPaPQYODOuAXH113Imccy5nsgJ++pW0DPijmd2SsX0O8KCZXVrBca2BF4EZwLBo8+nADsDvorOXzGMGAAMAWrVq1XHkyJE1yrxkyRKaN29eo2MrZMaWN93EZmPH8r9zz2XeIYfk7UfXSd46kqSskKy8ScoKycqbpKxQu7ydO3eebmY7ZX3QzAp2A74ABmXZvpRQVCo67iZgJrBW2ra1gVnArVX93o4dO1pNTZkypcbHVmrZMrPu3c0aNDAbPTpvP7bO8taBJGU1S1beJGU1S1beJGU1q11eYJpV8L5a6DaP98lo24i64TYloy0kw1bAf81sWWqDmf0E/Bf4RR3krHuNGsHIkbDjjnDkkfDmm3Encs65ait08XgGOFDSumnb+gLfA1MrOW4WsK2ktVMbJK0DbEs4I0mm5s1h7Fho1Qq6d4ePP447kXPOVUuhi8edwI/AE5L2i9olBgM3WVr3XUkfSRqedty9wGbAk5K6S+oBPAW0Bu4uWPq60KpVGAOyYkVYSGrhwrgTOedclQpaPMxsMbAv0BAYQxggeDMwKGPXRtE+qeOmAwcRBgo+BDxIuNS1v5klf86PX/86LGU7ezYccgh8/33ciZxzrlIFnxjRzGYAXarYpyzLtsnA5DqKFb899oARI+CII6BfP3jsMWhY5ZRfzjkXC18Mqpj07g033QRPPgnnneeDCJ1zRcunZC8255wTLl/dfHNYSOq88+JO5Jxza/DiUYxuuCGMQj///DCpYp8+cSdyzrnVePEoRg0awEMPhYWkjjkGWreGvfaKO5Vzzv3M2zyKVePG8PTTsMUWcOih8N57cSdyzrmfefEoZhtuGMaArL12GAOyYEHciZxzDvDiUfy22CKMQv/yyzAKfcmSuBM555wXj0TYaacw7uOdd0Lj+fLlcSdyzpU4Lx5J0a0bDBsWLmOddpqPAXHOxcp7WyXJySfDrFkwZEhYSOqyy+JO5JwrUV48kubKK8MYkD//ORSQ446LO5FzrgR58UgaCe65B+bNg5NOgs02g/33jzuVc67EeJtHEq29Njz+OGy9dZgP653kTyzsnEsWLx5J1aIFjB8fvnbrFi5lOedcgXjxSLI2bUIBWbIEdt8d2ralU5cuUFYWpnd3zrk64sUj6bbbDk4/HebOhTlzkFnokTVggBcQ51yd8eJRHzz88Jrbysth4MDCZ3HOlQQvHvXB7Nm5bXfOuVry4lEftGuX23bnnKslLx71wZAh0LTpmtt79y58FudcSfDiUR/06wd33w3t22MStG0bZuMdNgzeeCPudM65esiLR33Rrx/MnMnU558PbR2vvgqbbgo9esDHH8edzjlXz3jxqK9atYJnn4WVK8NCUgsXxp3IOVePePGoz7bcEkaPDqPPDz44dN91zrk88OJR3/3ud2Gw4Ouvh0tbK1bEncg5Vw948SgFvXrBLbfAU0/Buef6QlLOuVrzKdlLxVlnhWlLbroJ2reH88+PO5FzLsG8eJSS668P7R8XXBAmVezbN+5EzrmE8uJRSho0gAcfhPnz4dhjoXVr2HvvuFM55xLI2zxKTePG8PTT0KEDHHoozJgRdyLnXAIVvHhI2kbSZEnlkuZJukJSw2oe20vSm5K+l7RI0rOSmtV15npngw3gmWdgnXXCGJD58+NO5JxLmIIWD0ktgUmAAYcCVwDnA5dX49iTgIeBZ4CuwEnAh/ilt5opK4Nx42DRIujeHb77Lu5EzrkEKfQb76lAE6CXmX0LTJS0HjBY0tBo2xokbQTcDJxpZvekPfRknSeuzzp2hMceCwMI+/QJAwrXWivuVM65BCj0ZauuwISMIjGSUFA6VXJcn+jrA3UVrGR17Qp33hmmMjn1VB8D4pyrlkIXj62A99M3mNlsoDx6rCK7Ah8AJ0qaI2mZpNcl/a7uopaQk06Cyy6D++6DK6+MO41zLgEKXTxaAl9n2b44eqwimwK/Bv4EXAQcDCwFnpXUKt8hS9Lll8Nxx8GgQXD//XGncc4VOVkBL1NIWgb80cxuydg+B3jQzC6t4LjngP2Brmb2bLRtPWAWcLuZXZblmAHAAIBWrVp1HDlyZI0yL1myhObNm9fo2DjUJq+WLWO7Sy5h/X/9i3evuYbFO++c53SrK6XnttCSlBWSlTdJWaF2eTt37jzdzHbK+qCZFewGfAEMyrJ9KaGoVHTcI8BKoHHG9knAqKp+b8eOHa2mpkyZUuNj41DrvN98Y7b99mbNm5u9/XZeMlWk5J7bAkpSVrNk5U1SVrPa5QWmWQXvq4W+bPU+GW0bktoCTcloC8nwHqDottrhhKLi8mW99WD8eFh/fejWLSws5ZxzGQpdPJ4BDpS0btq2vsD3wNRKjhsbfe2c2iCpBdAReCffIUve5puHQYTl5aE31uLFcSdyzhWZQhePO4EfgSck7Re1SwwGbrK07ruSPpI0PHXfzKYBTwPDJR0nqTswGlgG/LWQf0DJ2HZbePJJ+PBD6NkTfvwx7kTOuSJS0OJhZouBfYGGwBjCyPKbgUEZuzaK9kl3NPAUcBPwOKFwdIl+pqsLnTuHnldTp0L//mFJW+ecI4apPcxsBtClin3KsmxbApwW3VyhHHVUaPe45BJo1w6uuy7uRM65IuDzQrmqXXRRKCBDh4YCcvrpcSdyzsXMi4ermgS33gpz5oQVCdu0CdO5O+dKlq/n4aqnUSP4xz/CZIq//z28/nrciZxzMaqyeEh6XtJW0ffHStqw7mO5otSsGYwdG1Yg7NEDPvoo7kTOuZhU58xjL2D96Pu/Ab+ouziu6G2ySZiB1yyMAfnyy7gTOediUJ02j8+AIyQtIYzo3iL6PquoN5Wrz371KxgzBrp0CWuBPP88NG0adyrnXAFVp3hcA9wBnENYAfDhCvZT9Hi1lpR1Cbf77vDww9C7d+jOO2oUNPR/eudKRZXFw8zukTQa+BXwInA64GcXLow8/8tfQg+ss8+G224LPbOcc/VelcVD0rHAODN7SdLlwNNmNq/uo7lEOPNMmDULbrwR2reHP/4x7kTOuQKoToN5eiP5n4E2dRfHJdLQoWEN9AsvhBqum+KcS5bqtHksBjaLvk+1azi3SoMG8MADsGBBWI2wdWvoVNmS9M65pKtO8ZgEPCTpg+j+/ZKWVrSzme2Sl2QuWRo3DrPw7rEHHHYYvPwybLNN3Kmcc3WkOsXjBMJkhFsBOwKfAt65361pgw3COiC77x7GgLz6Kmy2WdXHOecSpzq9rcqBGwEk7QcMNDNfgMllV1YG48bB3ntD9+7w4ouw7rpVHuacS5ac5rYysy28cLgq7bgjPP44vPsuHH44LFsWdyLnXJ7lPDGipA6Shkl6V9Lc6OsdkjrURUCXUAcdBHfdBc89B6ecEqYzcc7VGzlNyS6pIzAF+IGwrvjnQCugN9BPUmczeyvvKV0ynXhiWAfkiivCGJBBmQtGOueSKtf1PG4A3ga6Rm0hAEhqCoyPHq90lUBXYgYPDgVk8OCwkNTxx8edyDmXB7kWj12APumFA0KjuqQbgEfylszVDxLcfTfMnQsnnxx6Xx14YNypnHO1lGubx/dARet5bEC4nOXc6tZaKzSgb7ttaEB/++24EznnainX4jEOuFbSnukbo/vXAGPyFczVM+utB+PHQ8uW0LkztGlDpy5dQtfeESPiTuecy1GuxeM84BNgqqT5kt6RNB+YShg8eH6+A7p6ZLPN4PTT4ZtvYO5cZBYmVRwwwAuIcwmTU5uHmS0C9pR0ELAz0BqYD7xuZs/VQT5X3wwbtua28nIYOBD69St8HudcjeTaVbehma0ws2eBZ+sok6vPZs/ObbtzrijletlqrqShkraukzSu/mvXLrftzrmilGvxuBM4HPiPpNclDZC0Xh3kcvXVkCHZ1zvfZ5+CR3HO1Vyuc1sNNrMOwP7AB8BNwHxJI6JJE52rXL9+YdxH+/aYFM44dtwRHnwQnnoq7nTOuWrKeW4rADN73syOBTYFzgR+DUyQNFPSYEk+D7erWL9+MHMmU59/PvS2+uc/YZdd4Pe/h9deizudc64aalQ80uwE7E1Y62Mx8E/gJOAjSUfX8me7UtG0KYwZA5tvDgcfDB99FHci51wVajKrbntJgyR9DEwmdNc9AdjMzI4B2gN3AdfnNamr3zbeOCwkBWFG3i99vTHnillOxUPSFOBjQrF4COhgZgea2aNm9hOAma0AHibMtutc9f3qV+EMZO5c6NEjjP9wzhWlXM88vgC6AWVR4/msCvb7F7BFtgckbSNpsqRySfMkXSGpYXUDSGogaZokk9Qjx/yu2O22G/zjH/Dmm3DUUbBiRdyJnHNZ5Nrbqq+ZPWdW+co+ZrYsW2GR1BKYBBhwKHAFYUqTy3OIcRLQJof9XdIcdhjceis8/TScdZYvJOVcEcp1SnYAJLUBtgQaZz5mZuMrOfRUoAnQy8y+BSZG40QGSxoabavs97YEhgAXA/fWJLtLiDPOCKPOr78+LCR14YVxJ3LOpcl1epJ1gUeBA1Kboq/pHw0ruwTVFZiQUSRGAtcBnah6Vt4rgZcJDfWuvrv2WvjsM7joImjbNnTldc4VhVzbPK4B2gF7EQpHT2AfYDhhVt3dqjh+K+D99A1mNhsojx6rkKTtCQ31F+SY2SVVgwZw//2w997Qvz+88ELMgZxzKaqi+WL1naVPgD8RVgxcBuxqZm9Gj90ItDWzPpUcvwz4o5ndkrF9DvCgmV1aybFTCbP3XiipjFCsDjazsRXsPwAYANCqVauOI0eOrPbfmW7JkiU0b968RsfGIUl5q5u10XffscOZZ7LOwoW8ddttlG+RtS9GnauPz22xSFLeJGWF2uXt3LnzdDPbKeuDZlbtG7AU2Cv6/jvgwLTH9gW+ruL4ZcA5WbbPAa6u5LgjgQXAetH9MsKlsh7Vyd2xY0erqSlTptT42DgkKW9OWWfONNt0U7O2bc3mzq2zTJWpt89tEUhS3iRlNatdXmCaVfC+mutlq8+AjaLvPwTSu8ruStXL0C4GWmTZ3jJ6bA2S1iIMOLwOaCBpfSA1GWOzqB3G1Xft24eVCBcvhm7d4NtK+1Y45+pYrsVjIpCaAPFm4HRJr0SDB68EHqzi+PfJaNuQ1BZoSkZbSJpmhK65NxEKzGLgneixkYAviF0qdtghrIX+n/+EtdCXLYs7kXMlq1rFQ1ITSb0JgwRfkdTKzB4CehPaHhYDZxC60FbmGeDAjLOFvsD3hKVss1kCdM64pbrdXAr48nOl5MAD4Z57YOLEsHytjwFxLhZVdtWV1IEwsK8sbfO3kvqY2ZPAkzn8vjuBs4AnJF0HdAAGAzdZWvddSR8BU83sRDNbDryQkSmV5V0zez2H3+/qg+OPD2NABg8Ol7MGD447kXMlpzpnHkOBlYTuuU2B3xAuFd2V6y8zs8WEhvWGhDEdlxMufw3K2LURlY8XcaXuz3+GE06Ayy+H4cPjTuNcyanOIMHdgfPN7OXo/nuSTom+tjaz+bn8QjObAXSpYp+yKh6fyaoBiq4USXDnnWESxVNOCdO5H3RQ3KmcKxnVOfNoDXySse1jwpv3pnlP5Fx1rbUWPPYYbLcdHHEEvPVW3ImcKxnV7W3lrZKuOK27bujCu8EG0L07zJwZdyLnSkJ1i8cESV+kbkDqUtXk9O3RY84VVuvWYSGpH34IY0AWZx0y5JzLo+q0eeQyXbpz8dhmG3jqKTjggDCl+4QJ0HiNSZ+dc3lSZfEwMy8eLhk6dYIHHgiz7x53XFhUqkHOKy0756qhRut5OFe0jjwyTON+4YXQrl1YD8Q5l3dePFz9c8EFYRDhDTeEAnLmmXEncq7e8eLh6h8JbrkF5syBs8+GNm2gZ8+4UzlXr/gFYVc/NWwII0bArrvCUUfBq6/Gnci5esWLh6u/mjaF0aPDmcfBB8P//hd3IufqDS8ern7beOMwBkSCrl3hCx+K5Fw+ePFw9d8vfwljx8L8+eEMZOnSuBM5l3hePFxp2HVXGDkSpk0L40CWL487kXOJ5sXDlY5DDoHbboMxY+Css3whKedqwbvqutLyhz/ArFkwdGhYSOqii+JO5FwiefFwpeeaa8Io9IsvhrZtQ1de51xOvHi40tOgAfztb6EBvX//MCtv585xp3IuUbzNw5WmddaBJ5+ELbcMo8//85+4EzmXKF48XOlaf/2wkFSzZmEdkLlz407kXGJ48XClrV07GDcuLCDVrRt8+23ciZxLBC8ezv3f/8GoUTBjBhx+OCxbFnci54qeFw/nIKxAeM89MHEinHyyjwFxrgre28q5lP79wzoggwaFy1lXXBF3IueKlhcP59JddlkoIFdeGQrISSfFnci5ouTFw7l0EgwbFnpenXoqbL55mI3XObcab/NwLtNaa8Gjj8L228MRR8BVV0FZGZ26dIGysrDIlHMlzouHc9msu27owtu4cbiUNWsWMgvzYg0Y4AXElTwvHs5VpHVrWHvtNbeXl8PAgYXP41wR8eLhXGUWLMi+ffbswuZwrsh48XCuMu3a5bbduRJR8OIhaRtJkyWVS5on6QpJDas4ZmdJf5P0UXTcB5IGSWpcqNyuRA0ZAk2brr6tUaOw3bkSVtCuupJaApOAGcChwC+AGwlF7E+VHNo32vc64ENge+DK6GvvOozsSl2/fuHrwIHY7NmoeXP47jtYtCjeXM7FrNDjPE4FmgC9zOxbYKKk9YDBkoZG27K51swWpt1/QdIPwF2S2pvZrDrO7UpZv37Qrx9TX3iBffbaK3TfPeccaNMGevWKO51zsSj0ZauuwISMIjGSUFA6VXRQRuFIeTv6uln+4jlXhYYNQzfd3XYLReWVV+JO5FwsCl08tgLeT99gZrOB8uixXOwOrAQ+zk8056qpSRMYPTosYXvIIfC//8WdyLmCkxVw9lBJy4A/mtktGdvnAA+a2aXV/DmbAv8GxptZ/wr2GQAMAGjVqlXHkSNH1ijzkiVLaN68eY2OjUOS8iYpK6yZt/Hcuex4xhmsaNqUt26/nWUtW8aYbnVJf26LWZKyQu3ydu7cebqZ7ZT1QTMr2A1YBpyTZfsc4Opq/oy1gReBT4CW1TmmY8eOVlNTpkyp8bFxSFLeJGU1qyDv66+bNWlitvPOZkuWFDxTRerFc1ukkpTVrHZ5gWlWwftqoS9bLQZaZNneMnqsUpIEPAj8BuhmZlUe41yd2mUXeOQRmD4djjwSli+PO5FzBVHo4vE+GW0bktoCTcloC6nALYQuvoeaWXX2d67uHXww/PWvMHYsnHmmLyTlSkKhu+o+A/xR0rpm9l20rS/wPTC1sgMlXQKcAfQxs5fqNqZzOTr11DBp4rXXQvv2cPHFcSdyrk4V+szjTuBH4AlJ+0WN2oOBmyyt+240knx42v2jgKsJl6zmStot7bZxYf8E5yowZAgcdRRcconPuuvqvYKeeZjZYkn7ArcDY4CvgZsJBSQzV/qUJQdEX/tHt3THA/fnN6lzNdCgAdx3H8yfD8cfH2bl7dIl7lTO1YmCryRoZjOASv9HmVlZxv3+rFk0nCs+66wDTzwBe+4JPXvCSy/BdtvFncq5vPNZdZ3Lt/XXh2eegebNoVs3mDMn7kTO5Z0XD+fqQtu2MH48fPMNdO8O31Y0bZtzyeTFw7m68tvfwqhRMGMG9O4NP/0UdyLn8saLh3N1af/94d57YdIkOPlkHwPi6o2CN5g7V3KOOw4++wwuuyysQHjllXEncq7WvHg4VwgDB4ZBhFddFQrIySfHnci5WvHi4VwhSDBsGMydC6edBptvHnpiOawwInsAABo+SURBVJdQ3ubhXKE0agSPPhoa0vv0CZMpOpdQXjycK6TmzWHcONh449CF99NP407kXI148XCu0DbdNAwi/Okn6NoVFi2KO5FzOfPi4VwcttoqLGU7cyYceij88EPciZzLiRcP5+Ky557w0EPwyitwzDGwcmXciZyrNi8ezsXpiCPghhvg8cfhggviTuNctXlXXefidu65MHs23HxzWEjq7LPjTuRclbx4OBc3CW68MYxCP/dcaNMmzIXlXBHzy1bOFYOGDeHvf4fdd4d+/eDll+NO5FylvHg4VyyaNIGnnw6Xrg45BD74IO5EzlXIi4dzxWSjjcIYkEaNwhiQzz+PO5FzWXnxcK7YdOgAY8eGwtGjByxdGnci59bgxcO5YrTzzvDII/DWW9C3LyxfHnci51bjxcO5YtWjB9xxR5gL6/TTfSEpV1S8q65zxeyUU8I6INdcExrSL7007kTOAV48nCt+Q4aEMSADB4aFpI4+Ou5EznnxcK7oSTB8OMybByecAK1bw777xp3KlThv83AuCdZeG554An79a+jVC959N+5ErsR58XAuKVq0gPHjYd11wxiQOXPiTuRKmBcP55KkbdtQQL77LqyB/s03cSdyJcqLh3NJs/324RLWe++FCRR/+inuRK4EefFwLon23Rfuuw8mT4aTTvIxIK7gvLeVc0l1zDFhHZA//Sl04b3qqrgTuRJS8DMPSdtImiypXNI8SVdIaliN41pI+pukxZK+kTRC0oaFyOxc0br0Ujj55DAW5O67407jSkhBzzwktQQmATOAQ4FfADcSitifqjj8UWBL4CRgJXAd8BSwV13lda7oSWEKk7lz4bTTYPPNoXv3uFO5ElDoM49TgSZALzObaGZ3ApcD50lar6KDJO0OHAAcZ2ajzOxJ4GhgT0n7FSK4c0WrUaMwieIOO0DPntC6NZ26dIGyMhgxIu50lRsxAsrKkpE3SVmhzvMWunh0BSaY2bdp20YSCkqnKo773MxeTG0wszeAT6PHnCttzZuH0efLl8OCBcgszIk1YEDxvsmNGBHyzZpV/HmTlBUKkrfQDeZbAc+nbzCz2ZLKo8fGVHLc+1m2vxc95pwbOnTNXlfl5XDssXDeeeESF+T+tSbHVOfrRx+tOdV8eTn07w9XX13ln1tQ//tfcrJCxXkHDgzLHOdBoYtHS+DrLNsXR4/V5LgO2Q6QNAAYANCqVSteeOGFnIKmLFmypMbHxiFJeZOUFYo/b6fZs1GW7bZyJfN22w0gfAqFVUUmo9hU9nhtjs32+Mbvv5897/LlfLnxxlkeic/GM2YkJitUknf2bKbm6zVsZgW7AcuAc7JsnwNcXclxE4Gnsmz/O/BKVb+3Y8eOVlNTpkyp8bFxSFLeJGU1S0De9u3Nwlvz6rf27eNOll2S8iYpq1ne8gLTrIL31UK3eSwGWmTZ3jJ6LN/HOVc6hgyBpk1X39a0adhejJKUN0lZoSB5C1083iejjUJSW6Ap2ds0KjwuUlFbiHOlp1+/MNajfXtMCotH3X133q5x512S8iYpKxQkb6GLxzPAgZLWTdvWF/gemFrFcZtK2jO1QdJOhPaOZ+oiqHOJ1K8fzJzJ1Oefh5kzi/fNLSVJeZOUFeo8b6GLx53Aj8ATkvaLGrUHAzdZWvddSR9JGp66b2avAs8BD0rqJekwYATwkplNKuhf4JxzrrDFw8wWA/sCDQndci8HbgYGZezaKNonXV/C2cl9wIPAdKBnXeZ1zjmXXcEnRjSzGUCXKvYpy7Lta+D46Oaccy5GPiW7c865nHnxcM45lzNZCSwiI+lLYFYND98IWJjHOHUtSXmTlBWSlTdJWSFZeZOUFWqXt72ZZR1CXxLFozYkTTOzneLOUV1JypukrJCsvEnKCsnKm6SsUHd5/bKVc865nHnxcM45lzMvHlVL2tqeScqbpKyQrLxJygrJypukrFBHeb3NwznnXM78zMM551zOvHg455zLmRcP55xzOfPi4ZxzLmdePJxzzuWs4LPquvyIVmDsBgh4zMwWSWoDXAD8ApgJ3G1m78aXEiRdBIyPO0d1SWoCNDKz79K2bQycAWwDrAT+BdxhZt/Ek9K5+HlX3YgkEdYH6Q5sDWxAeKNYALwG3G9m/4sv4SqSdgEmAs2A5cBXwIHAeGAF8F9gW2BTYD8z+2dMUZG0EjDCcsEPAyPN7OO48lRF0njgQzM7O7q/O2G1ypWENWQEdAR+ArqY2X9jzLoD0MTMXknbdhBwCasK3TvA4PR9ikX0f+5gYEfCa2Qa4YNGUb8pSVqPMFdUFzN7Ke488HOmLsDawDgzWxp96DmdsOLqJ4QPk/Py9juL/N+pIKIneTzhTWEB4Y1hc8IL+hnCk/9r4EozuzKunCmSJhLOGnsCSwkLah1GeHM73MyWSVoHeApobGadY8y6ErgO2A7Yn5D7LUIhedTM5saVLRtJC4ETzezp6P5rhOf4sNTZiKQWwGjgBzM7MMasrwFjzGxIdP8E4F5gCvA8odDtC+wF9E79TTFlfYXwvL4X3W9JWB20I7Ak2q054YPagelnfnGQ9IdKHm4CXA/8BfgQwMzuKESubCT9EpgMtI02fQocQPiAuT7wMeH963ugo5nNycsvNrOSvwH/ILwItkvbthnwLDAqut+J8CI/oQjyLgK6pt3fhPAp84CM/boDC2POuhLYJfp+A2BA9EJfHt1eiLZtEPfzGmUsB/ZOu/9T5vOa9twujTnrt+nZgI+A27LsdyfwTrG8DqL7wwlnzAelbTsIWAzcXASvg5WEs/iVFdzSH1sRc9ZHCWeYv4z+jz0UvZ+9Aqwb7bNRtM9d+fq93mAedAUutrTr8hZO704FDpPU2symAlcDZ8eUMZNl+T7zNLKoTivN7Cszu9vM9gXaAOcTTrPvBOZLGhdrwOA/QPqZ2ueE/5CZNiQUmjitzLjfHng8y36PEz55FpNDgCvM7NnUhuj7IUCv2FKt8jTwBXAi0NDMGqRuhNeDgH2ibZlLZhfansAQM/vIzL4C/kRo97zBojM4M1sI3MLqr+1a8eIRNCB8ksi0gvAiaRHdfx3YslChKjENuEBSc0kNgEuBucBpkhoCSGoE/IHwZlh0zGyBmf3FzH4HbEFYx36zmGMBXAtcLOmE6DkcAlwvaX9Ja0taJ2pXuIbwiS9O/wT6pd3/L5Bt6u2dCa+PYrI+4XWcaTqhrS5WZtYTOA74I/CmpD3SH44nVYVaEi63p6T+rTPXMPqE8KEtL7y3VTARuErSv83sE/j5muythH+UVEN5c6AYetgMJGReTLj0U05oLHsM+FBSqsF8M8KlgKJmZrMIb9rXFkGWJySdSfiUdjPwAeHDw7OENw1Fu44mvLHE6VLg5egDxG2EhvIHJG1AuBwowuviHODiuEKm6S0pVdwWA9kWGdqIcDkudmb2nKTtCc/fOEnPEnozxtoek8UXhLPOlBXAXYSz5nSbkMfs3mAORF1cJxDOKmYRrnNvAfwI/N7Mnon2G0pYWatvXFlTosw9CB8ARpnZfEmbAhcSLlHMAu41s7dijImkQcA9lsdeHoUgaUOgL7AL4ZOwCG947wFjzWx6jPF+Jun/gGHArqxe3FLfLyZcHvpLPAmDqONEpvvN7ISM/e4CtjGzvQqTrHqi/1tDCZfU7iIUlM5m9mKswQBJTwFfZT6XWfa7DdjKzPbPy+/14hFEl3v6AL8FGhMaHx+OriE6V9QkbU0oIJmF7hUzWxZntlxIOhn42MyejztLNlHX7ZsJH9C6WxF0gZbUCmhqZp9Wsd95hI4Tk/Pye7141D+SGphZtk96RUNSY0Kj3krgo2J8g4vaPDqQNubHzGbHm8q54uAN5hkk/UZSb0knRbfekn4Td65MknpJekrSeEkHR9v6SpoJLJM0K/oUFytJR0fjD1L3G0m6ltBN89+EBv2vJBXDNXkAJHWUNJpwffg94GXgVeBTSXMlXSGpaawh6xFF4s6RjaQmmf/Wkv4vel/oGFeuohBn/+RiugEnENoJsvXtXkGY7uP4uHNGWftEuV4idCksB04mtNUMJ4wq/UeU+8CYs84ATku7f2OU9zJgD0I3w8GEAUyXFsFzewChrWsaoWv2YMJA0Z+izOcTejX9C2hZBHl7EMbNvBe9FvbOss+uxD8W4QCiMQdp2w4jDBhdDiyLnvPucT+nUbYWwJNRruXAPUBD4IGM94WXgI3izlvNv6l3Pl8Hsf9BxXADzoxeJH8ljMbdKHqhNIy+3xO4PXpTOb0I8r4J3Jl2v1+U7caM/f4GTIo5aznQKe3+F8DZWfa7AJhVBM/tdOCBCl4jMwln642jN707Ys66f/QG9nL0+pwe3b+R6JJ0tF8xFI8VrD5IsGf0BvxK9G9/PuHsbjlZBmXGkPdWwhQkZwLHRh8YRgGfRYVwY8L4sLnAsLjzVvNvymvx8DYPQNInhDfjoVXsdyFwqpl1KEyyCnN8C/Qys0nR/RaEBtL9LK2hMbqcdZeZxTZ+QtJ84AwzGxXd/5FwNvRCxn77A6PNrEnhU66W43vgEDObmLG9JWFk/2/M7D1JxwLXmVnrOHJGmV4izMN1fNq2EwhvfBMJPQV/kLQroeE8tsFsUW+r3czsjej+W8BcMzs4Y7/xQDMz6xRDzPQcMwkD7+6J7u9AKM7Hm9kDafudTDhj3iKWoCHDfdXctT1hYGNeXgfe5hG0Bt6oxn5vUAQDmAjdMNNfAKm5gb7O2G8JYTBWnEYTBjSuHd2fBPw+y36/J3y6i9sXhB53mX5LeN5T43xmsWrwaFy2Bf6evsHM7iNMpbMb8Hw05qMYbUvo8prpbsJEiXHbmFXjuyCaw4owT1S6j8g+XqWQjiOcDW1Xxa19RT+gJnyQYPAOcLKkF62CXkpRg97JhEbeuM0izJo6AcDMVkRdCN/L2O8XwPwCZ8t0CWEk9H8k3QuMAa6TtC2rBrJ1BnYgzLAat7uBKyU1JxS6nwgjtAcCU2zVeJUOQNw9r34gzKy8GjObHo2InkC4LDS4wLkqkn6Z4xuyD1hbSnF8qP2UUISnRvf3Ilxm+x2hnSNlD+J/HXwIvGFmx1a2k6TDgUfy9Uu9eATnE0YQz5D0BGH68NSn+BbAVoRrtG0ojhHbT5AxzYCZvZ5lvyNZ/YVecGb2laTdCG++5xFGuQLsHt1+Ilxi2cvM3own5SpmNiS6xHIx8OfUZkIHhHPSdl1GaFCP078J191HZz5gZp9EBWQ8cH+Bc1VkgqTl0fctCB8YpmbssxXxf+CBMN/aXyRtRyh0fQgfhC6LPli8QzhDOheIe6bt1whFrSrpg0hrzds8IpJ+QRidfRCrpjZO+YzQ4+Z6K+K1KDJJagd8bWZFMd0DgKQyVh/I9rEV5xiPtQhnbo2BT4rpOUyRdAphipIdrILBrJKaEXoN7WdhUr9YRDMNZPrQzB7O2O+FaHsxdDM/i3A5dS3CbA13Svo9oU0pNTHm3cBFcb6Goy7De5jZrVXstxGhzS6zYNfs93rxWFPUrzvVVvC1mcU9e6pzrkhEl7A3MrMv484SJy8e9Ux0Sv0W0K8YLgMpgcu6KiFL/DoXJy8eaaI3jc0Il1IWZnl8I6CbmT1Y8HCr5+hWycPNCI1iFxNNx25m4wuRKxslaFlXSNYSv9UVzXt1hJldEXOOgi+Vmk/RGUf6srnTCX9H7G+iCrMV9yb8f7rfzN6X9FvgclZ94LndzCbk7ZfGPXClGG7AOoTpzFdEt2WEkdotMvaLfbBVlCNJq5wtBA5Nu/8aYUT0umnbWhAaTicUwXM7kbCM6/qEa923A3MIo7fXSnu9PEPofRX767caf1NeB4fVMMMvCb0EU6/Ljwlvap8QCvSbhKnYPwfaFMFz9gqwddr9llHGlVHOb1k1yHHduHJG2Q4kfPhaED2v3xJ6MC6O8v01+n+3grCccn5+b9z/SMVwI/Sq+ZrQFXcn4KzoRfwh8Ku0/YqleEwj9Eg5ntB3O/22ffSi7pPaFnPWxCzrGuVI0hK/7ap5OzXu1y0xLZVai7yJWTaXMMPAY4QVDyF0olgMDM/Y7yHgtbz93rj/kYrhRuiae0bGtk2BF4Evgd2jbcVSPERY9/sLwpQJW6Q91iJ64a8xx1FMWd8ABqXd/ww4Mst+xwJfFkHeRRlvEBtHz+f+Gft1K4LikTrLrOpWDGeg84A+affbR7l6Zex3PPC/IngdZBaPL4FzsuwX+7Q6hK7E+6Xdbxnl75Kx3wGEDkB5+b0+ziNoS8bgPzNbIGlfQrWeJKkfxdH/HAuvhLslPQpcBbwbLfRyVbzJsroWGCHpM+BBVi3ruohwqSo1SLAYlnWFVUv8vkQ4a0pf4vd5CwMyi2WJ3++A54F7q9hvT0I39DjFslRqHhXzsrnfs/pg0dT3mVP9NCUMLM0LLx7BPOBXhDONn1nou32kpL8QTgtjbSjPZGZfA2dIupvQ9/xD4DqKaI1lS9ayrpCsJX7fILTLjatsp2jtlLjFslRqLSVl2dyXgT9L+jDKcgNhNuuLolkzvovmv7uQUOzywntb8fPEYh3MbJ9K9rmE8KnZLMYJ5ioj6UjCUpltCBOgxb5EZooSsqwrJGqJ38uAAWaWOag1c7+9gcvNrHNhkmXNEMtSqTWlBC2bK+mXhKl0Uq+DmYSz+ccJI/ZnAWWED0Odzexfefm9Xjx+7ubWF7jGKll2VtJRhGvfx1e0T9yiSyrNgCVmtiLuPM5BfEul1jUVybK50fiuPQg9BCeb2ffRYOeTWPWB52Ezm5O33+nFwznnXK6KYfZKV0ck3SNpeNw5qiNJWSF5eZ3LN28wz4Gke4AGZnZi3FmqqTPJ+YCQpKyQoLySJhGuMuwbd5aqJCkrJCtvvrN68chNYt4wAMzsl3FnqK4kZYXE5RXJed0mKSskK29es3qbRz0WddHcxMziXqymSknKCsnL61y+JaViFgVJjaM1MpKiO2FFtCRIUlZIUF5JayXldZukrJCsvPnO6sUjN4l5w3ClQdLpkj6W9L2kdyQdk2W3HSmC122SskKy8saR1ds8EkhSdfuUZxsRW1BJygrJyhsNCr2NsETu24SlSO+XdChwtJnlbSqK2kpSVkhW3riyevEgWW8Ykb0J03zMqGK/YpiWIklZIVl5LwBuMLOf562K5mMbAUyR1MPMFsWWbnVJygrJyhtLVm8wByQtp3pvGJsDu8Y9PYmkd4D3zaxvFfsdDjwSZ94kZY1yJCavpO+Ag83shYztZYT1RhoS5t/aGHjFs1ZfkvLGldXbPIL/Av8xsyMquwE3xR008hqwWzX2S594MC5JygrJyvsNYWK+1ZjZTMKli4XAq8DOhY2VVZKyQrLyxpLVzzz4eXKzg8ysfRX79SasaR1r0ZX0C+A3Zja6iv2aELqTZk57XTBJyhrlSExeSU8D35nZ0RU83oQwOV5XYp7QM0lZozyJyRtXVi8eJOsNw7kUSUcA5wI9KprQU1JDYBhhQs8tCpkvI0diskZZEpM3rqxePJxzzuXM2zycc87lzIuHc865nHnxcCVHUn9J0yV9J2mxpLcl1UlPOklbShosaf1q7DtYkqXd5kkaFbXJVXXs/ZKyrbHtXJ3w4uFKisJywvcCE4BewLHA08AhdfQrtwQGAVUWj8g3wO7R7QLg/4DJkppVcdyVQP8aZnQuZz7C3JWaM4C7zOzStG1jJF0eV6AMy83stej71yTNBv4JdAMey9xZUhMz+97MPi5kSOf8zMOVmvWBBZkbLa3boaSy6LLRUZIeii5vfSFpUOZxkrpIel3SD5I+l3SHwnrSSNoHGBPt+mn0M2fmmHd69LUs+pkzJd0o6TJJc4Bvo+1rXLaS1F7SPyQtlFQu6d+Sjkp7vLGkoZI+k/RjNKFetxzzuRLlZx6u1LwFnBl9oh9bxZw/1wNjgcMJc14NkrTQzP4KIOk3wLPARKA30Ba4FuhAmA7iLaJ5hwiXyOYDP+aYtyz6ml7wjiLMivAHKvg/LGkTwqji8ijDZ8C2UcaUx4FdCJfVPgb6AKMl7WRm/8oxpysxXjxcqTkdeAq4HzBJ7wGjCBPLfZux73/N7JTo+wnRG/KlkoaZ2UrgMmAWcIiZrQCQ9BXwiKTdzexVSR9Ex78dTRdRJUmp/5cdgDuA74BJGbv1qGK21HOBFkBHM5sfbZuc9jv2JSwxsI+ZTY02PydpS2AgcER1srrS5ZetXEkxs38DWxMayO8gzE91GTAtdbkpzZMZ958ANgPaRPd3AZ5MFY7IKGA5sGcNI24ILItuHxAKSN+0AgAwuRrTbHcBns04Lt1+hLOZlyU1St0IBWanGmZ3JcTPPFzJMbMfCW0RYwAknUjogXUi8Je0Xb/IODR1vzUwO/r6ecbPXiFpEbBBDeN9Q3hjN8Kb+zxbcxqIz9c4ak0bAm9W8vhGwKaEIpVpRZZtzq3Gi4creWY2XNJQYKuMhzap4P78tK+r7RPNIbQhkHWOoWpYbmZVjdeozpxCiwjFrSJfAXOBw6obzLl0ftnKlZSo3SJz28aE9oHMT/Q9M+6nGr3nRPdfB3pGBSN9n0bAS9H9n6KvhV48ajJwoKRWlTy+KbDEzKZl3goX0yWVn3m4UvNuNIX1c4TLUO0JvZHKgQcy9v1NNF3/KEJvqxOBs6PGcoCrCMt+PiVpGKEt5Dpggpm9Gu2TajA/RdJIoNzM3q2bP201NxMGQP5T0hBCb6utgWZmNpTQQ2wCMFHSdYTeW+sRBiU2NrNLCpDRJZgXD1dqrgAOBW4ltEssAF4hNEp/mrHvhUAPQvH4gTCK+/bUg2b2X0ldgasJjenfEtaRvjBtn1mSLgDOAs4knLWU1cUfls7MvpS0BzAUuAVYB/gQuCZ63CT1Ai4FzgHaES5l/YuwHrZzlfIp2Z3LEC3f+Slhac+x8aZxrjh5m4dzzrmcefFwzjmXM79s5ZxzLmd+5uGccy5nXjycc87lzIuHc865nHnxcM45lzMvHs4553L2/z15xUlOCG0UAAAAAElFTkSuQmCC\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",
    "# construct amplitude estimation \n",
    "ae = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha, \n",
    "                                  state_preparation=european_put,\n",
    "                                  objective_qubits=[num_uncertainty_qubits],\n",
    "                                  post_processing=european_put_objective.post_processing)"
   ]
  },
  {
   "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.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)"
   ]
  },
  {
   "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.1801\n",
      "Confidence interval:\t[0.1742, 0.1860]\n"
     ]
    }
   ],
   "source": [
    "conf_int = np.array(result['confidence_interval'])\n",
    "print('Exact value:        \\t%.4f' % exact_value)\n",
    "print('Estimated value:    \\t%.4f' % (result['estimation']))\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 factory 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",
    "# construct amplitude estimation \n",
    "ae_delta = IterativeAmplitudeEstimation(epsilon=epsilon, alpha=alpha, \n",
    "                                        state_preparation=european_put_delta,\n",
    "                                        objective_qubits=[num_uncertainty_qubits])"
   ]
  },
  {
   "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.run(quantum_instance=Aer.get_backend('qasm_simulator'), shots=100)"
   ]
  },
  {
   "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",
      "Esimated value: \t-0.8180\n",
      "Confidence interval: \t[-0.8244, -0.8117]\n"
     ]
    }
   ],
   "source": [
    "conf_int = -np.array(result_delta['confidence_interval'])[::-1]\n",
    "print('Exact delta:    \\t%.4f' % exact_delta)\n",
    "print('Esimated 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>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>None</td></tr><tr><td>Terra</td><td>0.16.0.dev0+28d8c6a</td></tr><tr><td>Aer</td><td>0.6.1</td></tr><tr><td>Ignis</td><td>0.5.0.dev0+470d8cc</td></tr><tr><td>Aqua</td><td>0.8.0.dev0+ce81016</td></tr><tr><td>IBM Q Provider</td><td>0.8.0</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.7.7 (default, May  6 2020, 04:59:01) \n",
       "[Clang 4.0.1 (tags/RELEASE_401/final)]</td></tr><tr><td>OS</td><td>Darwin</td></tr><tr><td>CPUs</td><td>2</td></tr><tr><td>Memory (Gb)</td><td>16.0</td></tr><tr><td colspan='2'>Fri Oct 16 10:45:35 2020 CEST</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, 2020.</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 qiskit.tools.jupyter\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.7.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
  }
 },
 "nbformat": 4,
 "nbformat_minor": 1
}