quantum_unary_option_pricing.ipynb

{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Quantum unary approach to option pricing\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _Abstract_\n",
    "\n",
    "Here is presented a novel quantum algorithm for European option pricing in finance, where the key idea is to work in the unary representation. The algorithm is divided in two main parts, amplitude distribution using a low depth circuit and computation of the expected return computed with simple controlled gates. Unary representation offers several advantages and one disadvantage. On the positive side, unary representation simplifies the structure and depth of the quantum circuit and also provides a post-selection consistency check that allows for a substantial mitigation in the error of the computation. On the negative side, unary representation offers lower resolution given a fixed number of qubits, as compared to binary algorithms.\n",
    "\n",
    "In the following the unary algorithm is described and simulations are perfomed with the inclusion of a simple noise model in order to test the robustness of the algorithm."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Introduction\n",
    "\n",
    "The use of quantum computation to solve financial problems, such as pricing of financial derivatives, is being discussed more and more nowadays [[1](https://arxiv.org/abs/1807.03890)], [[2](https://arxiv.org/abs/1805.00109)], [[3](https://arxiv.org/abs/1905.02666)]. Algorithms that estimate the expected payoff of a derivative have been recently presented, and the expectation is that, provided a robust quantum device, they will outclass current classical means. However, for Noisy Intermediate Scale Quantum (NISQ) devices currently available, these algorithms are too costly to run at high scales.\n",
    "\n",
    "One of the most promising advantages provided by quantum computation is the exponential Hilbert space at its disposal. In the case of option pricing this directly translates to high precision in its estimations. However, the circuits that load the expected probabilities, and more importantly, the comparator needed to gauge if the option price is higher or lower than the accorded strike quickly becomes prohivitive in NISQ devices.\n",
    "\n",
    "Here, a different point of view is proposed. Changing to the unary representation yields a simpler circuit and allows a post-selection scheme that results in error mitication. The unary representation means states where only a single qubit is in the $|1\\rangle$ state. This represents a loss of Hilbert space of $2^n\\longrightarrow n$, but is beneficial in other fronts. This allows for a simple amplitude distribution algorithm to encode the probability distriution of the asset price into the quantum register, and a more straight forward way to compute the expected payoff."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Pricing of European options\n",
    "\n",
    "European options, known also as path-independent _vanilla_ options are the some of the easiest to approach at first. This is due to the fact that they only depend on the asset price at maturity time $T$, refrenced henceforth as $S_T$. The expected value of said asset at time $T$ can be evaluated via the Black-Scholes model [[4](http://www.jstor.org/stable/1831029)]. The Black-Scholes model assumes that the price of the asset underlying the option is dictated by a geometric Brownian motion. The evolution of the stock price can be described by the stochastic differential equation\n",
    "\n",
    "$$dS_t=S_t r dt+S_t\\sigma dW_t,$$\n",
    "\n",
    "where $r$ and $\\sigma$ are the interest rate and volatility of the asset and $W_t$ characterizes a Brownian motion. This Brownian motion can also be written as $dW_t=\\sqrt{dt}N(0,1)$ with $N(0,1)$ a normal distribution with $0$ mean and standard deviation of $1$. This stochastic equation can be solved exatly for constant interest rate and volatility as\n",
    "\n",
    "$$S_t=S_0e^{\\sigma W_t+(\\alpha-\\sigma^2/2)t}.$$\n",
    "\n",
    "This is an exponential of a normal random variable, hidden in $W_t$, therfore the probability distribution function of $S_t$ will take on a _log-normal_ form. At a maturity time T the probability distribution function of the asset price becomes\n",
    "\n",
    "$$p(S_T)=\\frac{1}{S_T\\sigma\\sqrt{2\\pi T}}\\exp{\\left(-\\frac{(\\ln S_T-(\\alpha-\\frac{\\sigma^2}{2})T-\\ln S_0)^2}{2\\sigma^2T}\\right)},$$\n",
    "\n",
    "from which $(\\alpha-\\frac{\\sigma^2}{2})T+\\ln S_0)$ is usually given as $\\mu$.\n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Lognormal probability distribution function normalized for representation in finite intervals\n",
    "\n",
    "'''\n",
    "def log_normal(x, mu, sig, T):\n",
    "    dx = x[1]-x[0]\n",
    "    log_norm = 1 / (x * sig * np.sqrt(2 * np.pi * T)) * np.exp(- np.power(np.log(x) - mu, 2.) / (2 * np.power(sig, 2.) * T))\n",
    "    f = log_norm*dx/(np.sum(log_norm * dx))\n",
    "    return f"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Once the probability distribution is computed the payoff of an European option can be defined as \n",
    "\n",
    "$$\\text{payoff}(S_T)=\\max(0,S_T-K)$$\n",
    "\n",
    "for call options, and \n",
    "\n",
    "$$\\text{payoff}(S_T)=\\max(0,K-S_T)$$\n",
    "\n",
    "for put options, where $K$ is the accorded strike price when the option was first aquired. An option provides the chance to buy ( _call_ ) or sell ( _put_ ) the underlying asset at the strike price $K$ at maturity time $T$. As it is not mandatory to do so, the option holder will not buy or sell the asset if there is no profit to be made."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Difficulties for NISQ devices in the binary representation\n",
    "\n",
    "Once the probability distribution function is succesfully uploaded into the binary representation, assumed to be efficient due to the integrability of the Black-Scholes model [[5](https://arxiv.org/abs/quant-ph/0208112)], in order to diferentiate states with $S_T>K$ ( _for call options_ ) a quantum comparator has to be implemented. This comparator [[6](https://arxiv.org/abs/quant-ph/0410184)], tailored to the strike price $K$, computes the carry bit in the most significant qubit to determine if the option price is above or below the strike. This procedure is based on multiple CNOT and Toffoli gates. Toffoli gates are a problem for NISQ devices as this three-quibit gates have to be decomposed to multiple single and two-qubit gates \n",
    "\n",
    "<img src=\"Images/Toffoli.png\" width=\"600\" title=\"Toffoli decomposition\"/>\n",
    "\n",
    "as seen in [7]. This makes the comparator hard to implement in near-term devices as the gates needed to complete the algorithm has a large decomposition in basis gates.\n",
    "\n",
    "The unitary representation does not need of such a comparator, therefore no three qubit gate will be needed."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Algorithm in the unary representation\n",
    "\n",
    "#### Advantadges for NISQ devices\n",
    "\n",
    "Unlike the binary case, there is no need of Toffoli gates for the comparator in the unary representation. At most, the most complicated gates that appear are CNOT and $cR_y$ gates. Overall, the probablity loading algorithm scales with $(n-1)$ partial-SWAP gates, where $n$ is the qubit number. If this number were to be decomposed into regular CNOT and single qubit gates, we would find $4(n-1)$ CNOT gates and $(2n-1)$ single qubit gates.\n",
    "\n",
    "The payoff circuit is even simpler. Instead of relying on a comparator and a payoff calculator, in the unary representation they can be applied at the same time with $(n-k)$ $cR_y$ gates, where $k$ refers to the qubit number where the option price becomes higher than the strike. Therefore, the number of $cR_y$ gates needed is always lower than the number of qubits, sometimes significantly so.\n",
    "\n",
    "The exact decomposition of the partial-SWAP gates and $cR_y$ gates will be shown in the following sections."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Amplitude distribution algorithm\n",
    "\n",
    "This algorithm distributes the mplitude of one qubit with the rest in order to represent a target probability distribution. The circuit is initialized with only the middle qubit occupied. We can introduce partial-SWAP gates\n",
    "<img src=\"Images/SWAPRy2.png\" width=\"200\" title=\"SWAP-Ry\"/>\n",
    "which can be decomposed as\n",
    "<img src=\"Images/SWAPRy.png\" width=\"250\" title=\"SWAP-Ry decomposition\"/>\n",
    "using only CNOT and $cR_y$ gates. In turn, in order to represent it in the native gates of the device, this $cR_y$ gate can be separated into\n",
    "<img src=\"Images/CRyDef.png\" width=\"300\" title=\"C-Ry decomposition\"/>\n",
    "\n",
    "The purpose of this two qubit gates is to share the probability of a state with its neighbours. The angle in which the rotation is performed encodes the amount of probability that leaves the occupied state to go to the next one. As this parameters can be trained or tuned analytically, one only needs one of this gates for each two qubit pair of the system, amounting to $(n-1)$ partial-SWAP gates. The completed amplitude distributor circuit for $n$ qubits is\n",
    "<img src=\"Images/Prob_text.png\" width=\"600\" title=\"Random walk circuit\"/>\n",
    "where the gates are applied starting from the centre, as it is the only occupied state when the algorithm starts.\n",
    "\n",
    "After this circuit is applied one finds the desired distribution encoded into the amplitudes of the system as\n",
    "$$|\\Psi\\rangle=\\sum_i^n\\sqrt{p_i}|q_i\\rangle.$$"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Circuit that emulates the probability sharing between neighbors seen usually\n",
    "in Brownian motion and stochastic systems\n",
    "\n",
    "'''\n",
    "def rw_circuit(qubits, parameters):\n",
    "    C = QuantumCircuit(qubits+1)\n",
    "    if qubits%2==0:\n",
    "        mid1 = int(qubits/2)\n",
    "        mid0 = int(mid1-1)\n",
    "        C.x(mid1)\n",
    "        #partial-SWAP gates\n",
    "        #-------------------------------------------------\n",
    "        C.cx(mid1, mid0)\n",
    "        C.cu3(parameters[mid0], 0, 0, mid0, mid1)\n",
    "        C.cx(mid1, mid0)\n",
    "        #-------------------------------------------------\n",
    "        for i in range(mid0):\n",
    "            #partial-SWAP gates\n",
    "            #-------------------------------------------------\n",
    "            C.cx(mid0-i, mid0-i-1)\n",
    "            C.cu3(parameters[mid0-i-1], 0, 0, mid0-i-1, mid0-i)\n",
    "            C.cx(mid0-i, mid0-i-1)\n",
    "            #-------------------------------------------------\n",
    "            C.cx(mid1+i, mid1+i+1)\n",
    "            C.cu3(parameters[mid1+i], 0, 0, mid1+i+1, mid1+i)\n",
    "            C.cx(mid1+i, mid1+i+1)\n",
    "            #-------------------------------------------------\n",
    "    else:\n",
    "        mid = int((qubits-1)/2) #The random walk starts from the middle qubit\n",
    "        C.x(mid)\n",
    "        for i in range(mid):\n",
    "            #partial-SWAP gates\n",
    "            #-------------------------------------------------\n",
    "            C.cx(mid-i, mid-i-1)\n",
    "            C.cu3(parameters[mid-i-1], 0, 0, mid-i-1, mid-i)\n",
    "            C.cx(mid-i, mid-i-1)\n",
    "            #-------------------------------------------------\n",
    "            C.cx(mid+i, mid+i+1)\n",
    "            C.cu3(parameters[mid+i], 0, 0, mid+i+1, mid+i)\n",
    "            C.cx(mid+i, mid+i+1)\n",
    "            #-------------------------------------------------\n",
    "    return C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "##### Parameter optimization\n",
    "\n",
    "Each parameter of the partial-SWAP gates for the random walk circuit can be solved exactly if the probability distribution function is analytic, as is the case in this example.\n",
    "\n",
    "Due to the step-like nature of the circuit one can define a system of $n$ equations with $n-1$ free parameters to be fulfilled simultaneously. Then, starting from the edges of the circuit, one can start solving the equations. Once the center qubit is reached the rest of necessary solutions can be computed directly, as shown in the following code."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Solving for the exact angles for the random walk circuit that enables the loading of\n",
    "a desired probability distribution function\n",
    "\n",
    "'''\n",
    "def rw_parameters(qubits, pdf):\n",
    "    if qubits%2==0:\n",
    "        mid = int((qubits)/2)\n",
    "    else:\n",
    "        mid = int((qubits-1)/2) #Important to keep track of the centre\n",
    "    last = 1\n",
    "    parameters = []\n",
    "    for i in range(mid-1):\n",
    "        angle = 2 * np.arctan(np.sqrt(pdf[i]/(pdf[i+1] * last)))\n",
    "        parameters.append(angle)\n",
    "        last = (np.cos(angle/2))**2 #The last solution is needed to solve the next one\n",
    "    angle = 2 * np.arcsin(np.sqrt(pdf[mid-1]/last))\n",
    "    parameters.append(angle)\n",
    "    last = (np.cos(angle/2))**2\n",
    "    for i in range(mid, qubits-1):\n",
    "        angle = 2 * np.arccos(np.sqrt(pdf[i]/last))\n",
    "        parameters.append(angle)\n",
    "        last *= (np.sin(angle/2))**2\n",
    "    return parameters"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Circuit to sample the created probability distribution funtion\n",
    "\n",
    "'''\n",
    "def measure_probability(qubits):\n",
    "    C = QuantumCircuit(qubits+1, qubits)\n",
    "    C.barrier(qubits)\n",
    "    C.measure(range(qubits),range(qubits)) #No measure on the ancilla qubit is necessary\n",
    "    return C"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "From the retuned sampling, extract only probabilies of unary states\n",
    "\n",
    "'''\n",
    "def extract_probability(qubits, counts, samples):\n",
    "    form = '{0:0%sb}' % str(qu)\n",
    "    prob = []\n",
    "    for i in range(qu):\n",
    "        prob.append(counts.get(form.format(2**i), 0))\n",
    "    prob = np.array(prob)\n",
    "    prob = prob/sum(prob)\n",
    "    return prob"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 6,
   "metadata": {},
   "outputs": [],
   "source": [
    "def extract_probability_all_samples(qubits, counts, samples):\n",
    "    \"\"\"\n",
    "    From the retuned sampling, extract only probabilies of all states\n",
    "    \"\"\"\n",
    "    form = '{0:0%sb}' % str(qubits)\n",
    "    prob = []\n",
    "    for i in range(2**qubits):\n",
    "        prob.append(counts.get(form.format(i), 0))\n",
    "    prob = np.array(prob)\n",
    "    prob = prob/sum(prob)\n",
    "    return prob"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Simulation\n",
    "\n",
    "In order to test that this first step of the algorithm is correct, we can run a simulation using *qasm_simulator*"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 7,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Import auxiliary packages\n",
    "\n",
    "'''\n",
    "import matplotlib.pyplot as plt\n",
    "%matplotlib inline\n",
    "import numpy as np\n",
    "from math import pi"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 8,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Import Qiskit packages\n",
    "\n",
    "'''\n",
    "from qiskit import QuantumCircuit, execute\n",
    "from qiskit import Aer"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 9,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Get qasm backend for the simulation of the circuit\n",
    "\n",
    "'''\n",
    "backend_sim = Aer.get_backend('qasm_simulator')"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 10,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Define all the parameters to be used in the computation\n",
    "\n",
    "'''\n",
    "qu = 8 #Number of qubits for the probability distribution (ancilla not included)\n",
    "S0 = 2\n",
    "K = 1.9\n",
    "sig = 0.4\n",
    "r = 0.05\n",
    "T = 0.25\n",
    "mu = (r-0.5*sig**2)*T+np.log(S0) #mu value to be fed to lognormal distribution\n",
    "samples = 10000"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to sucesfully plot our data, it is useful to get the mean and variance of the log-normal probability distribution. The distance of 3 standard deviations from the mean is enough to represent the relevant part of the distribution."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 11,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Set the relevant zone of study and create the mapping between qubit and option price, and\n",
    "generate the target lognormal distribution within the interval\n",
    "\n",
    "'''\n",
    "mean = np.exp(mu+0.5*T*sig**2)\n",
    "variance = (np.exp(T*sig**2)-1)*np.exp(2*mu+T*sig**2)\n",
    "\n",
    "S = np.linspace(max(mean - 3 * np.sqrt(variance), 0), mean+3*np.sqrt(variance), qu)\n",
    "\n",
    "ln = log_normal(S, mu, sig, T) "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 12,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Solve for the parameters needed to create the target lognormal distribution\n",
    "\n",
    "'''\n",
    "lognormal_parameters = rw_parameters(qu, ln)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 13,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Build the probaility loading circuit with the adjusted parameters\n",
    "\n",
    "'''\n",
    "prob_loading = rw_circuit(qu, lognormal_parameters)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 14,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Circuit to test the precision of the probability loading algorithm\n",
    "\n",
    "'''\n",
    "circ_prob = prob_loading+measure_probability(qu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 15,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Run the test circuit through the simulator and sample the results\n",
    "in order to get the estimated probabilities\n",
    "\n",
    "'''\n",
    "job_sim = execute(circ_prob, backend_sim, shots=samples)\n",
    "result_sim = job_sim.result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 16,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_sim = result_sim.get_counts(circ_prob)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 17,
   "metadata": {},
   "outputs": [],
   "source": [
    "prob_sim = extract_probability(qu, counts_sim, samples)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 18,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Generate a very precise target probability distribution to benchmark the quantum results\n",
    "\n",
    "'''\n",
    "Sp = np.linspace(max(mean - 3 * np.sqrt(variance), 0), mean+3*np.sqrt(variance), 10000)\n",
    "lnp = log_normal(Sp, mu, sig, T)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 19,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "image/png": "iVBORw0KGgoAAAANSUhEUgAAAagAAAEYCAYAAAAJeGK1AAAABHNCSVQICAgIfAhkiAAAAAlwSFlzAAALEgAACxIB0t1+/AAAADh0RVh0U29mdHdhcmUAbWF0cGxvdGxpYiB2ZXJzaW9uMy4xLjEsIGh0dHA6Ly9tYXRwbG90bGliLm9yZy8QZhcZAAAgAElEQVR4nOzdd3wUdfrA8c+TnhB6kyqItABSDEVEFFFAUcB24s8GttPDdt7p6XmWs9zZG6iInl0PEKWoFLGgIh0pGpAqvdcA6cnz+2Mm3BJSNmQ3s9k879drXzs79ZmZnX12vvOd74iqYowxxoSaCK8DMMYYYwpjCcoYY0xIsgRljDEmJFmCMsYYE5IsQRljjAlJlqCMMcaEJEtQ5UREmorIYRGJDIFYrhaRr8phOSoip7rdo0XkoQDN95htKSKzROSmQMzbnd80Ebk+UPPzmW+8iHwuIgdF5JNAzz+YRGSYiMwuZnhQtlm4EJG/i8hbQZp3QL//ocQSVBHcA/IXEUkTkR0i8rqI1CjF9BtE5Lz8z6q6SVUTVTU3OBH7T1U/UtV+5bzMW1X18ZLGK7jdiphXwLaliDwqIh8WmP8FqvpeWeddiMuB+kBtVb0iEDMUkT+IyEoROSQiK0RkSCDmW1q+26ykZFYZqeq/VNXzJOLP8RVKLEEVQkT+AjwN3AtUB3oAJwMzRSTGy9jKSkSivI6hLCp4/CcDq1U1p7QTFrbeItII+BC4B6iG8339WETqlTVQYwry5NhTVXv5vHAO9MPAHwr0TwR2Aze4nx8FJgDjgEPAz0BHd9gHQB6Q7s7rPqAZoECUO05DYAqwD1gL3OyzrEeB8cD77rxTgORiYlbgTmA9sAd4Fohwhw0DfgJeBPYCT7j9ZvtM3w6Y6cayE/i72z8CuB9Y5047HqhVTBz3AtuBbcANblynusPeBZ5wu+sAXwAH3GX+6C6ruO12I7AJ+KGQbTkL+DewAEgFJufHCZwDbCkQ5wbgPGAAkAVku8tb5jO/m3y2wT+AjcAud59Ud4flx3G9G9se4MEits0/CyzrRj/nfXS9C5lnd2BXgX67gTOKiKE2zncu1d1Wj+d/Dwpu00K2wzCc79Eo4CDwG9C34LhAWyADyHXX84A7/EJgBc73eSvw1zIco9e622wv8GD+/nSHdQPm4ny3trvxxrjDBOc42OVug1+A9j7fz9eAaW7cPwEnAS8B+9317exHbH9z1+8QsCp/G+Ec0x8W2NbDgc3u/G8FugLL3dhHFfg9+NDn8zH7qsB+agF8626bPcBHQI2ifpfc/oNwfmMOuPNqW+BY+ZsbV6bv96Ncfo/Lc2EV4YXzo5VT2I4A3gP+6/OlycYptokG/gr8DkT77NjzivlS/eAeEHFAJ5wflnN95p2Bc1BH4vz4zismZgW+A2oBTYHVHPvDkgPcAUQB8fgkKKAqzoH8FzeWqkB3d9hdwDygMRALvJG//kVst51Ae6AK8DFFJ6h/A6Pd7RYNnAVICdvtfXe+8YVsy1k4Pwr5y/6U//0YnEMRCcpnW39YYPgsn+13A84fiFNw/qR8BnxQILY33bg64hzEbYvYRscsy895H13vQuYXCXyP8wMTCQwBtgBVilj+WJw/GVXcbbWV0iWoHODP7j67EidR1Spi3NkFlr0dOMvtrgl0OcHjMwnnx7U3znfyBTeu/P15Ok6JR5S7TiuBu91h/YHFQA2cZNUWaODz/dzjTh+H8yP/O3Cdu22fAL4rIbbWOAmnoc82bVFw3/ts69HusvrhHO+TgHpAI5wkenYR35tj9lWBbX8qcL67beri/M68VNh33/3cCjjiThON86dwLf9L6huApUATCvkOBvtlRXzHqwPs0cKLYba7w/MtVtUJqpqNc6DE4RwcxRKRJsCZwN9UNUNVlwJv4RwM+War6lR1rrN8gPPjV5ynVXWfqm7C+dd3lc+wbao6UlVzVDW9wHQXATtU9Xk3lkOqOt8ddivOGcEWVc3EOVAuL+JU/w/AO6r6q6oeccctSjbQADhZVbNV9Ud1j4ZiPKqqRwqJP98HPst+CPhDgCqkXA28oKrrVfUw8AAwtMA2+KeqpqvqMmAZJe+r0sy7yPV2vxvv4/wZyHTf/+hug2O42+Iy4GF3fr/i/OEqjV04P3bZqjoO5wxhoJ/TZgNJIlJNVfer6s+lXHa+y4EvVPUH9zv5EM5ZAQCqulhV57nf9Q04f6rO9omhKtAG5w/RSlXd7jPvie70GcBEIENV33e38zigcwmx5eIkhiQRiVbVDaq6rpjxH3ePua9wksR/VXWXqm7FKVUoaXnHUdW1qjpTVTNVdTfO79LZxUxyJfClO0028BzOn62ePuO8oqqbizn2gsYS1PH2AHWK+BFu4A7Ptzm/Q1XzcP69NvRjGQ2Bfap6yKffRpx/Tvl2+HSnAXEllAFv9uneWCCOzRStCU4RXmFOBiaKyAEROYDzbzQX50J/QQ0LiaEoz+L8S/tKRNaLyP3FjJuvuHUoOHwjzr/BOkWMWxoNOXZdNuL8O/fdBgX3VWIA513kersXu5/BOUuMwfkhektEOhUyel133v7uo8JsLfBHouD3rDiX4ZQIbBSR70XkjMJGEpEUt4bmYRE5q5BRjvmeucl4r8/0rUTkC7diUyrwL9zvgap+i1Pk9yqwS0TGiEg1n3nv9OlOL+RzsftVVdcCd+P8OdslImNFpLjtU6blFUZE6rvL3equ/4cUfxwc8x10f8c2c+xvUUnHXtBYgjreXJx/o5f69hSRROAC4Buf3k18hkfgFIVtc3sVd0awDaglIlV9+jXFKXI5UU18upv6xFFSLJtxipiKGnaBqtbwecW5//AK2l5IDIVyz9L+oqqn4BRP3SMifUuItaQzrILLzsb5M3EESMgf4J5J1C3FfLfhJGrfeedw7I/JifJn3sXF1wnn2tQiVc1T1YXAfJzrawXtdudd1D7KP+tK8Ol3UoF5NBIRKTD9No53XMyqulBVB+MUYU3CKWo8fkLVdurU0ExU1R8LGeWY75mIJOBcW8v3Os71opaqWg34O05xXv78X1HV03GKClvhXDcNGFX9WFV74exXxalsVVbHfIc5fr/4+pe73A7u+l+Dz/pz/L455jvo7t8mHPtbVNIxEjSWoApQ1YM4F7RHisgAEYkWkWY4B9QWnOK2fKeLyKXumc3dOIltnjtsJ0X88KvqZmAO8G8RiROR03Auhn9Y2Ph+uldEarrFh3fhFEn44wuggYjcLSKxIlJVRLq7w0YDT4rIyQAiUldEBhcxn/HAMBFJcn80HilqgSJykYic6h4MB3HOyvKLaYrcbiW4xmfZjwET3KKZ1ThnnwNFJBqnUkKsz3Q7gWbuH4zC/Bf4s4g0d/+k/AsYV0QRcGmVdd4LgbPyz5hEpDPO9bzlBUd0t8VnwKMikiAiSTiVO/KH78b5UbpGRCJF5AacC+6+6gF3usfEFTjXcKYWEtdOoHF+jVcRiRHn3rvqbjFSKj7FcqU0AbhIRHq583+MY3/HqrrzPywibYDb8geISFcR6e5+D47gXPc50TiOIyKtReRcEYl1550eoPkvBXqLc/9fdZyi4KJUxblGd9Ct5VkwARc8vsYDA0Wkr7td/oLzOzYnAHGXmSWoQqjqMzj/vJ7D+bLPxzmb6OuWe+ebjFOGux+nZtGl7gEITkWAf7jFY38tZDFX4Vzs3IZT3v2Iqn5dhrAn41wAXgp8CfzHn4ncYsbzgYtxiqrWAH3cwS/j1Pr6SkQO4STf7kXMZxrOta9vcYrvvi1msS2Br3EOpLnAa6r6nTuspO1WlA9wLnTvwLkWeKcb10HgTzjX+Lbi/DBt8Zku/4bZvSJS2HWRt915/4Bz0TwDp8JJIJRp3qr6PW5tUnf/fAr8y72mUZjbcYqNduBsq3cKDL8Z5wdtL07NzoI/UvNx9t0e4EngclXdy/G+xakVtkNE8ovErwU2uMVOt+Jcfys1VU0BRuBcb9uOc+z57s+/Av+HU4vuTY79o1bN7bef/9UCfPZE4ihCLPAUzvbZgZPQi0smflHVmTjrsRznGP+imNH/CXTB+eP3Jc6fEl/HHF+qugrnLGukG/fFwMWqmlXWuAMhv+aUKSUReRSnhto1IRCL4hRprPU6FlNxiMgwnNpfvbyOpSxEZAPOepTlD54JQXYGZYwxJiRV5LvyjTGmXIlIU5wbjguT5N7mYQLEiviMMcaEJCviM8YYE5LCpoivTp062qxZM6/DMMYYU4LFixfvUdW6JY0XNgmqWbNmLFq0yOswjDHGlEBE/GrFxIr4jDHGhCRLUMYYY0KSJShjjDEhKWyuQRljTKBkZ2ezZcsWMjIyvA6lQouLi6Nx48ZER0ef0PSWoIwxpoAtW7ZQtWpVmjVrxrENuBt/qSp79+5ly5YtNG/e/ITmYUV8xhhTQEZGBrVr17bkVAYiQu3atct0FmoJyhhjCmHJqezKug0tQRljjAlJdg3KmEDLSoPfvoTfZ8GOXyBtP0REQOJJ0LAznNoXTukDkXb4maJFRkbSoUOHo5+HDh3K/fffH5B5L126lG3btnHhhRcGZH7BYkeIMYGSkQpzXoH5b0BmKsTXchJS3baguXBwCyx+F+a/7iSrM++E5BshOs7ryE0Iio+PZ+nSpUGZ99KlS1m0aFHIJygr4jMmENbMhJGnww/PQotzYdiXcO86uPYzuPQNuOwtuGE63L8RrvwQ6raGGX+HV7vC+lleR28qiIMHD9K6dWtWrVoFwFVXXcWbb74JwG233UZycjLt2rXjkUceOTrNwoUL6dmzJx07dqRbt24cPHiQhx9+mHHjxtGpUyfGjRtX6LJCgZ1BmUrv4pGzix64Z02x0wp5XMuXXCHfsEEb8DL3sDalKaSsAYqZts59nFatP7elvkrj9wczMe4S3k0YTp5E8vkdFfoBtyZA0tPT6dSp09HPDzzwAFdeeSWjRo1i2LBh3HXXXezfv5+bb74ZgCeffJJatWqRm5tL3759Wb58OW3atOHKK69k3LhxdO3aldTUVBISEnjsscdYtGgRo0aN8mr1/GIJypgTFEUOd/FfzpHFTNOevMklZOP/DYnLoztyV42RDD/yNpdkTOTk3I08k/i3IEZsTsg7A4Mz3+FfFju4qCK+888/n08++YQRI0awbNmyo/3Hjx/PmDFjyMnJYfv27axYsQIRoUGDBnTt2hWAatWqBXYdgiyoCUpEBgAvA5HAW6r6VIHhtwIjgFzgMHCLqq4QkWbASmCVO+o8Vb01mLEaU6g6LQvtHaG5/OXwM/TKWsx7CdczIe4KOIEqtVkSyxuJt7E+6hRuO/IaT6b+HdK6QUKtskZuwlReXh4rV64kISGB/fv307hxY37//Xeee+45Fi5cSM2aNRk2bFhYtIIRtAQlIpHAq8D5wBZgoYhMUVXfxyV/rKqj3fEHAS8AA9xh61S1E8aEGlX+dORVemXN5j8JNzIp/tIyz3JmXH/2RtTmwUNPwLsXwfVToEqdAARryqyEM53y9uKLL9K2bVv+9a9/MXz4cObOnUtqaipVqlShevXq7Ny5k2nTpnHOOefQunVrtm/fzsKFC+natSuHDh0iPj6eqlWrcujQIa9XpUTBrCTRDVirqutVNQsYCwz2HUFVU30+VgHs+fMm5A3JmEj/zBmMj/9DQJJTvp9jkvlntUdh3zr4+A+QdSRg8zYVT/41qPzX/fffz6pVq3jrrbd4/vnnOeuss+jduzdPPPEEHTt2pHPnzrRp04b/+7//48wzzwQgJiaGcePGcccdd9CxY0fOP/98MjIy6NOnDytWrKjUlSQaAZt9Pm8BuhccSURGAPcAMcC5PoOai8gSIBX4h6r+WMi0twC3ADRt2jRwkRtThNOylzEs7R3mxPTkg/jrAj7/5dGd4PK3Ydw1MOEGuPIju1+qksrNzS20/8qVK492v/DCC0e733333ULH79q1K/PmzTuu/8KFC8sWYDnwvJq5qr6qqi2AvwH/cHtvB5qqamec5PWxiBx3dU9Vx6hqsqom161b4tODjSmTGnn7uffQ02yNbMRLiX8+oWtOfmkzEC58FlZPh5kPBWcZxlQAwUxQW4EmPp8bu/2KMhYYAqCqmaq61+1eDKwDWgUpTmNKpsqIw6NI0DSeSvw76ZIQ3OV1vQm63wrzXoNfPw3usowJUcFMUAuBliLSXERigKHAFN8RRMS3itRA3BtHRKSuW8kCETkFaAmsD2KsxhSrT9a39MiexwcJ17E5qpyKk89/HJr0gMl3wK6VJY9vTJgJWoJS1RzgdmAGTpXx8aqaIiKPuTX2AG4XkRQRWYpTlHe92783sNztPwG4VVX3BStWY4pTI28/txx5g1+j2jElbnDJEwRKVAxc8S7EVHGuR2VX/GrDxpRGUK++qupUYGqBfg/7dN9VxHSfAlauYULC9WnvEquZjEq8kzznxL78VGsAQ16Hjy6Dbx+H/k+W7/KN8ZDnlSSMCWWts1dyXubXTI4bwtbIxt4E0fI855rU3Ffh9x+8icEYD1j9VWOKEKG53HbkdfZE1GZcwlBvgzn/MVj3HUy8Df40F+IqVpM1FV2x7TWeAH/bW9yyZQsjRoxgxYoV5ObmcuGFF/L8888TGxsbsFhmzZpFTEwMPXv2DNg8A8XOoIwpwtlZs2iRu453Em4gQ+K9DSamClzyBqRuhW+f8DYWUy5UlUsvvZQhQ4awZs0a1qxZQ3p6Ovfdd19AlzNr1izmzJkT0HkGiiUoYwoRpdlcnfYRayNb8GNMb6/DcTTp6hT1LRgDWxZ7HY0Jsm+//Za4uDiGDx8OOA8wfPHFF3n//fcZNWoUt99++9FxL7roImbNmgUU/diNZs2a8cgjj9ClSxc6dOjAb7/9xoYNGxg9ejQvvvginTp14scff2TYsGFMmDDh6HSJiYmAk8jOPvtsBg8ezCmnnML999/PRx99RLdu3ejQoQPr1q0L+DawBGVMIfpnTKd+3k7eT7gelRA6TPo+DFVPgs/vgtwcr6MxQZSSksLpp59+TL9q1arRrFkzcnKK3vdPPvkkixYtYvny5Xz//fcsX7786LA6derw888/c9ttt/Hcc8/RrFkzbr31Vv785z+zdOlSzjrrrGJjWrZsGaNHj2blypV88MEHrF69mgULFnDTTTcxcuTIsq1wIULoyDMmNMRqBlemj+WXqA4sie7idTjHiqsGFzwDO39xbuI1poDx48fTpUsXOnfuTEpKCitW/K997ksvddqOPP3009mwYUOp5921a1caNGhAbGwsLVq0oF+/fgB06NDhhOZXEktQxhRwYcaX1NQDvJdwffCaMyqLthdDqwHw/dNwaKfX0ZggSUpKYvHiY4tyU1NT2bFjB7Vr1yYvL+9o//xHa+Q/duObb75h+fLlDBw48JjHbuRXroiMjCzyLCwqKurovPPy8sjKyjpueoCIiIijnyMiIoo9qztRlqCM8RGtWQxJn8iS6E6sim7rdTiFE4H+/4KcTPj2Ma+jMUHSt29f0tLSeP/99wGn8di//OUv3H777TRv3pylS5eSl5fH5s2bWbBgAUChj90oScFHbzRr1uxoYpwyZQrZ2dlBWDv/WDVzY3z0zfyGWrqf5+Pv9TqU4tVuAd3/6Nwb1fVmaGiPTgsmf6uFB5KIMHHiREaMGMHjjz/O7t27ufLKK3nwwQdRVZo3b05SUhJt27alSxenKNr3sRtNmjQ5+tiN4lx88cVcfvnlTJ48mZEjR3LzzTczePBgOnbsyIABA6hSpUqwV7VIohoej2BKTk7WRYsWeR2GqYDy73GJ0FxGH7iF1Ihq/LXaC54V7/n9Y5hxEF7p4jz1d/i00CyOrKBWrlxJ27ahdQY9Z84crrrqKiZOnHg0IVUEhW1LEVmsqsklTWtFfMa4zsyaTYO8HSf8+PZyF1cd+j4Em+ZCykSvozFB1rNnTzZu3FihklNZWYIyBkCVy9I/ZXNkY+bH9PA6Gv91vhbqt4dvHoNc764VGBMMlqCMAZJyVtAidx2T44aE1n1PJYmIhL6PwP7f4ef3vY4mrITL5Q8vlXUbVqAj0ZjguSjjcw5LFWbF9vE6lNJreT407elUO89K8zqasBAXF8fevXstSZWBqrJ3717i4uJOeB5Wi89UerVy99Az6yemxA0mU078YPKMCJz3CLzdH+aPhrPu8TqiCq9x48Zs2bKF3bt3ex1KhRYXF0fjxif+FABLUKbSuyBzGoLyZdxAr0M5cU17QKsL4KeXIHk4xNf0OqIKLTo6mubNm3sdRqVnRXymcsvJpH/GdBZGd2VnZAOvoymbvg9BRirMftHrSIwJCEtQpnJLmURNPcAXcYO8jqTs6reDDlfAgjfhsBVNmYrPEpSp3H5+j20RDVgW3dHrSALj7PsgJwPmBr5laWPKW1ATlIgMEJFVIrJWRO4vZPitIvKLiCwVkdkikuQz7AF3ulUi0j+YcZpKau862PgTM+P6Vayq5cWp0xLaXwYL3oIje72OxpgyCdpRKSKRwKvABUAScJVvAnJ9rKodVLUT8AzwgjttEjAUaAcMAF5z52dM4Cz5ACSSb2P7eh1JYPW+F7LT7CzKVHjB/NvYDVirqutVNQsYCwz2HUFVU30+VgHybzoYDIxV1UxV/R1Y687PmMDIzYGlH0PLfuyLqO11NIFVtzW0vxTmj7GzKFOhBTNBNQI2+3ze4vY7hoiMEJF1OGdQd5ZmWmNO2Jqv4PBO6HKd15EER+/73LOoUV5HYswJ87zgXVVfVdUWwN+Af5RmWhG5RUQWicgiu6HOlMrP70NifWjZz+tIgqNeG2g3BBaMgbR9XkdjzAkJZoLaCjTx+dzY7VeUscCQ0kyrqmNUNVlVk+vWrVvGcE2lcWiHcwbV8SqIDON71XvfB1mHYd7rXkdizAkJZoJaCLQUkeYiEoNT6WGK7wgi0tLn40Bgjds9BRgqIrEi0hxoCSwIYqymMvnlE9Bc6HyN15EEV/0kaHORcxaVedjraIwptaAlKFXNAW4HZgArgfGqmiIij4lI/l2Rt4tIiogsBe4BrnenTQHGAyuA6cAIVc0NVqymklk+Dhq6D/oLd2feDRkH4Of3vI7EmFILavmGqk4Fphbo97BP913FTPsk8GTwojOV0q7fYMcvMOApryMpH026QrOzYM4o59HwUTFeR2SM3zyvJGFMufplPEikczNrZXHm3XBom7PuxlQglqBM5ZGX51x/OuUcSKzndTTl59S+UL8D/PSysw2MqSAsQZnKY/N8OLAJTvuD15GULxHodTfsWQ2rppY8vjEhIozr2BpTwC/jISoe2lTc5z5dPHL2CU0XoXUZHXESqZ/+k79Wq+4krTL6/I5eZZ6HMcWxMyhTOeRkQcpEJznFVvU6mnKXJ5F8Fn8ZrXNW0z7nF6/DMcYvlqBM5bDuG0jfX/mK93x8E3se+6UGl6d/4nUoxvjFivhM5ZAy0XkMeotzvY4kOPasKXGUbOBzzuS67C9puucHNuHHE4Qrw71iJmTZGZQJfzmZsGqaU7wXGe11NJ6aRk8yNIYhfOd1KMaUyM6gTPhb9y1kpkLSJV5HEjx+nukcBr45fD79Mmfwfs07OBBRK7hxGVMGdgZlwl/KJIirAaec7XUkIWFy/BAiyWVgxpdeh2JMsSxBmfBmxXvH2R7ZkPnRPbgwYyqxmuF1OMYUyRKUCW/rZ0HmQUgaUuKolcmk+CFU01TOzfzG61CMKZIlKBPeUiZBXHWneSNz1IqodqyObMXg9EmIWvNHJjRZgjLhKycLVn0JrQdaK94FiTApfgiN8rbRNdsetWZCkyUoE75+/x4yDjqPPjfH+SmmF7si6jIkfZLXoRhTKEtQJnylTILYala8V4Q8iWRK3GA65PxCi5ySb/Q1prxZgjLhKTcbfvsCWl8IUbFeRxOyZsb2I03iuSR9otehGHMcS1AmPG2a6zzqvO3FXkcS0tIiqjAjdgC9sn6kTu5ur8Mx5hiWoEx4WjUdImOhRR+vIwl5n8cNAuDijCkeR2LMsSxBmfCjCqunQfPeEFPF62hC3u7IevwU04v+mdOJ1zSvwzHmqKAmKBEZICKrRGStiNxfyPB7RGSFiCwXkW9E5GSfYbkistR92V874789a2DfemjV3+tIKoxJcUOoommcl/G116EYc1TQEpSIRAKvAhcAScBVIpJUYLQlQLKqngZMAJ7xGZauqp3c16BgxWnC0OrpznurAd7GUYGsiW7Niqi2DMqYRITmeh2OMUBwz6C6AWtVdb2qZgFjgcG+I6jqd6pHyxTmAY2DGI+pLFZPh/odoEYTryOpUCbHXcJJeTvpnjXP61CMAYKboBoBm30+b3H7FeVGYJrP5zgRWSQi80Sk0DstReQWd5xFu3dbDSQDpO2DTfOseO8EzIvpwY6I+gzOmOx1KMYAIVJJQkSuAZKBZ316n6yqycD/AS+JSIuC06nqGFVNVtXkunXrllO0JqSt/QY0F1pf4HUkFU6eRPJ53CDa5aTQMnuV1+EYE9QEtRXwLWNp7PY7hoicBzwIDFLVzPz+qrrVfV8PzAI6BzFWEy5WT4Mq9aBhF68jqZBmxvbjiCQwOMOaPzLeC2aCWgi0FJHmIhIDDAWOqY0nIp2BN3CS0y6f/jVFJNbtrgOcCawIYqwmHORmw5qvoVU/iAiJwoEKJz0iga9i+9Mra7bduGs8F7SjWFVzgNuBGcBKYLyqpojIYyKSXyvvWSAR+KRAdfK2wCIRWQZ8BzylqpagTPE2zXWe/WS198ok/8bdizI+9zgSU9lFBXPmqjoVmFqg38M+3ecVMd0coEMwYzNhaPUMiIyBU6z1iLLYHVmPOTFn0j9zOmMTriJD4r0OyVRSVg5iwscqt/WI2ESvI6nwJsUNIVGPcF7GTK9DMZWYJSgTHvasgX3rrHgvQFZHt2FlVBsGZUy2G3eNZyxBmfBwtPUIu/8pUCbFXUKDvB10y5rvdSimkrIEZcLDqulQvz3UaOp1JGFjXswZ7Iyob1XOjWcsQZmKL32/U4PPzp4Cyrlx92La56TQMme11+GYSsivBCUin4nIQBGxhGZCz5qvndYjWlnrEYH2VWx/0iSewel2FmXKn78J5zWcJofWiMhTItI6iDEZUzqrp27SdBcAACAASURBVENCHWh0uteRhJ30iARmxA7gTLtx13jArwSlql+r6tVAF2AD8LWIzBGR4SISHcwAjSlWbjasnekU71nrEUHxRdzFCMpAu3HXlDO/j2gRqQ0MA27CeY7TyzgJy26UMN7ZNA8yrPWIYNoVWZ+5MT0ZkDmdOE33OhxTifh7DWoi8COQAFysqoNUdZyq3oHTVJEx3lg93Wk9osW5XkcS1ibFXWI37ppy5+8Z1JuqmqSq/1bV7QD5jbm6j8Qwxhurp0Ozs6z1iCBbFd2G36LacHHGFLtx15QbfxPUE4X0mxvIQIwptT1rYe9aK94rJ5PihtAwbzvdshd4HYqpJIptLFZETsJ5Cm68+2gMcQdVwynuM8Y7q90HMLe2BFUe5sb0ZGdEPQanT2RezBleh2MqgZJaM++PUzGiMfCCT/9DwN+DFJMx/lk9A+q1s9Yjykn+E3dvSnuLU3PWAL28DsmEuWKL+FT1PVXtAwxT1T4+r0Gq+lk5xWjM8dL3w8Y51npEOZsZ289u3DXlpqQivmtU9UOgmYjcU3C4qr5QyGTGBN/ab5zWI1pb6xHlKS2iCjNj+zEw4ws4uBWqN/I6JBPGSqokUcV9TwSqFvIyxhvWeoRnPo8bhKCwYIzXoZgwV+wZlKq+4b7/s3zCMcYPuTmwZia0vhAiIr2OptLZGXkSc2POoNfid6D3vVbF3wRNSUV8rxQ3XFXvDGw4xvhh8zzIOGC19zw0OW4IvVJ/gqUfQ/dbvA7HhKmSavEtLpcojCmNVdOs9QiP/RbVFholw/zXoeuNdiZrgqKkIr73yjJzERmA02ZfJPCWqj5VYPg9OG375QC7gRtUdaM77HrgH+6oT5Q1FhNGVs+AZr0g1i6DekYEzhgBE4Y71wPbDPQ6IhOGiq0kISIvue+fi8iUgq8Spo0EXgUuAJKAq0QkqcBoS4BkVT0NmAA8405bC3gE6A50Ax4RkZqlXz0Tdvaug71rrPWIUNB2EFRvAnNf9ToSE6ZKKuL7wH1/7gTm3Q1Yq6rrAURkLDAYWJE/gqp+5zP+POAat7s/MFNV97nTzgQGAP89gThMOFnlth5hCcp7kVHQ/Y/w1T9g2xJo2NnriEyYKelG3cXu+/c4be/tB/YBc91+xWkEbPb5vMXtV5QbgWmlmVZEbhGRRSKyaPdue5hapbB6OtRLgponex2JAehyHcQkwtzXvI7EhCF/H7cxEFgHvAKMAtaKSMDukBSRa4Bk4NnSTKeqY1Q1WVWT69atG6hwTKhKPwCb5lrrEaEkrjp0vhZSPoPUbV5HY8KMv62ZPw/0UdVzVPVsoA/wYgnTbAWa+Hxu7PY7hoicBzwIDFLVzNJMayqZtV9DXg60stYjQkr3P4Lm2Y27JuD8TVCHVHWtz+f1OA3GFmch0FJEmotIDDAUOKZihdtC+hs4yWmXz6AZQD8RqelWjujn9jOV2erpkFAbGtsjyEJKrebQ5iJY9DZkHvY6GhNGSqrFd6mIXAosEpGpIjLMrf79OU4CKpKq5gC34ySWlcB4VU0RkcdEZJA72rM4zSh9IiJL82sGupUjHneXsRB4LL/ChKmk8luPaNnP7rkJRWeMgIyDsMzqMZnAKakW38U+3TuBs93u3UB8STNX1anA1AL9HvbpPq+Yad8G3i5pGaaS2DzfaT3Cau+FpibdnXYR570GyTfYnwgTECXdqDu8vAIxplirp0FEtLUeEaqO3rh7A/z2BSQN9joiEwZKOoMCQETicKqBtwPi8vur6g1BisuYY+W3HhFXzetITFHaDoaazeHHF5ybeEVKnsaYYvhbSeID4CScG2i/x6lVV1IlCWMCY+862LPaivdCXWQU9Lobti+F9d+VPL4xJfA3QZ2qqg8BR9w28QbiNENkTPCtnu68W+vloa/jVVC1gXMWZUwZ+Zugst33AyLSHqgO1AtOSMYUsGoa1G0LNZt5HYkpSVSscy1qw4+wZZHX0ZgKzt8ENca9H+khnHuZVgBPBy0qY/JZ6xEVz+nDIK6GnUWZMvOrkoSqvuV2fg+cErxwjClg3TdO6xGtrfWICiO2qtO6xPdPw66VUK+t1xGZCsrftvhqi8hIEflZRBaLyEsiUjvYwRnDqukQXwsad/U6ElMa3W+F6ASY/ZLXkZgKzK8zKGAs8ANwmfv5amAcUOSNtsYU6x0/HnCnClvmQ3xNeG9QyeMDDP+ybHGZwEio5RT1zX8D+vzdWp83J8Tfa1ANVPVxVf3dfT0B1A9mYMaQmeoU78XbyXqFdMbtIBEwZ6TXkZgKyt8zqK9EZCgw3v18OdZ4qynBxSNnFzP0gRKnH3bkbQaxkqtzHyf9cIJ/Cy1imZ/f0cu/6U3gVG8EHYfCkg/g7Psg0Sr+mtIpqbHYQyKSCtwMfAxkua+xwC3BD89UZl2zF/BrdHvSI/xMTib0nHk35GTC3FFeR2IqoJKeqFtVVau57xGqGuW+IlTV2pwxQdMgdxtNczezMLqb16GYsqhzKnS4HBa8BUf2eB2NqWD8vQaFiAwSkefc10XBDMqYblnzAVgQYw2WVHi974PsNLsWZUrN32rmTwF34dyguwK4S0T+HczATOXWLWs+GyJPZmfkSV6HYsqqbitofxkseBOO7PU6GlOB+HsGdSFwvqq+7T6naQBOe3zGBFyVvEO0y0mxs6dwcrZ7FjXXzqKM//wu4gNq+HRXD3QgxuRLzl5MJHksiLYEFTbqtob2l8L8MXYWZfzmb4L6N7BERN4VkfeAxcCTwQvLVGbdsuazX2qwOqqV16GYQMq/FmU1+oyfSkxQIiLAbKAH8BnwKXCGqo4LcmymEorSbLpkL2ZhTDdUSnOCb0JevTbQ7hJYYGdRxj8l/gKoqgJTVXW7qk5xXzvKITZTCSXlpJCoR1gQY9XLw9LZ90HWETuLMn7x9y/qzyJS6tY6RWSAiKwSkbUicn8hw3u7DdDmiMjlBYblishS9zWltMs2FVO3rAVkEsPS6M5eh2KCoV5bO4syfvM3QXUH5onIOhFZLiK/iMjy4iYQkUjgVeACIAm4SkSSCoy2CRiG00pFQemq2sl9+dlSqKnQVOmeNY9l0Z3IlDivozHBcvbfnLOo2fa8KFM8f9viO5GnxXUD1qrqegARGQsMxrmPCgBV3eAOyzuB+Zsw0zR3Iyfl7WRC/BVeh2KCqV4bp42+BW9Cjz85bfYZU4iS2uKLE5G7gXtx7n3aqqob818lzLsRsNnn8xa3n7/iRGSRiMwTkSFFxHeLO86i3bt3l2LWJhR1y14AwEK7/hT+znkANM95qKExRSipiO89IBn4Baeo7vmgR/Q/J6tqMvB/wEsi0qLgCKo6RlWTVTW5bt265RiaCYbuWfNZHdmKfRH2eI2wV/NkSB4OSz6EPWu9jsaEqJISVJKqXqOqb+A8YuOsUsx7K9DE53Njt59fVHWr+74emAXYVfMwViNvP61yVlntvcqk970QFQvf2S2VpnAlJajs/A5VzSnlvBcCLUWkuYjEAEMBv2rjiUhNEYl1u+sAZ+Jz7cqEn+SshUSgzLfmjSqPxHrQ4zZI+Qy2L/M6GhOCSkpQHUUk1X0dAk7L73afE1UkN6HdjvNgw5XAeFVNEZHHRGQQgIh0FZEtwBXAGyKS4k7eFlgkIsuA74CnVNUSVBjrnjWPXRF12RDZ3OtQTHnqeSfE1YBvHvc6EhOCiq3Fp6qRZZm5qk4Fphbo97BP90Kcor+C080BOpRl2abiiNUMOmcvYUbcABDxOhxTnuJrQK+74etHYeMcOLmn1xGZEGJtyRjPnZ61mFiymBtjP06VUrc/QtUG8NVDoOp1NCaEWIIynjsjaw4HpRorogrex20qhZgEOPcfsHUR/Pqp19GYEGIJyngqSrPpmr2A+TE9yJMylSibiqzjVVC/A3z9T8jO8DoaEyIsQRlPdcxeRhVNY44V71VuEZHQ/wk4uAnmv+51NCZEWIIynuqRNYc0iWdZdCevQzFeO+UcaDUAfnwBDlvLMMYSlPFQhObSI2seC6O7kiPRXodjQsH5jzsNyc76t9eRmBBgCcp4pm3OSmroQau9Z/6nbivoeiMsfhd2/eZ1NMZjlqCMZ87ImkMW0SyOSfY6FBNKzr4fYhLhqwet2nkl5+/jNowJLFXOyJrDkujOZEi819GY8vLOQP/GS6gFa7+G186ABD8aDx7+ZdniMiHJzqCMJ1rkrqVe3m4r3jOFq9oAohNg33rIy/U6GuMRO4MynuiZNYdcIlhgjcOGpYtHzi5iyAN+z6N9/HL+nfoA/93TjI8Tril+5CKXB5/f0cvvZZrQYmdQpvypclbmDyyP7sihiGpeR2NC1K/Rp/F9TG8uS5/ASbnbvQ7HeMASlCl3p+aupUHeDn6MKc3jxUxl9HbCjeQQxc1HxngdivGAJShT7s7K/IEcIu36kynRvsg6jE24im7ZC+iatcDrcEw5swRlypcqvbJ+ZEl0Zw5HVPU6GlMBfB43iE2RTfjjkdHEqrXTV5lYgjLlqnXOb9TL282Psb29DsVUEDkSzetVRlA/bydXp33odTimHFmCMuXqrKwfySaK+dE9vA7FVCC/RndgeuwABmVMpkXOGq/DMeXEEpQpN6J59MqczeLoZNIiqngdjqlg3k0YzgGpwR2HXyFSc7wOx5QDS1Cm3CTlrKC27uXHWKu9Z0rvSEQio6vcRovc9QzJmOh1OKYcWIIy5easzB/IJNZuzjUnbF5sT+bE9OSqtI9pkLvN63BMkAU1QYnIABFZJSJrReT+Qob3FpGfRSRHRC4vMOx6EVnjvq4PZpwm+CI1hzOzZrMwJtna3jNl8kaVW8mWaO44/AqieV6HY4IoaAlKRCKBV4ELgCTgKhFJKjDaJmAY8HGBaWsBjwDdgW7AIyJSM1ixmuDrnP0zNfQgs2L7eB2KqeD2RdTmrYSb6ZDzC4MyJnsdjgmiYJ5BdQPWqup6Vc0CxgKDfUdQ1Q2quhwo+DeoPzBTVfep6n5gJjAgiLGaIOuT+R2pUo3F0fZoDVN238Sex7zo7lyX9h5NcjZ5HY4JkmAmqEbAZp/PW9x+AZtWRG4RkUUismj3bntEdKhKyDtC96x5/BDT256cawJDhFcT7yBd4rnn8HNWqy9MVehKEqo6RlWTVTW5bt26XodjinBm1mxiyeI7K94zAXQgoiavVbmdU3PXcWX6WK/DMUEQzAS1FWji87mx2y/Y05oQ0yfzO7ZENGJ1VGuvQzFhZk7smXwbcy5/SB9Hy+xVXodjAiyYCWoh0FJEmotIDDAUmOLntDOAfiJS060c0c/tZyqYerk76ZDzC9/FngsiXodjwtCYKn9kX0Rt/nr4WeLz0rwOxwRQ0BKUquYAt+MklpXAeFVNEZHHRGQQgIh0FZEtwBXAGyKS4k67D3gcJ8ktBB5z+5kK5pzM7wCs9p4JmiMRiTyXeC/183Yy4shIUPU6JBMgQX2irqpOBaYW6PewT/dCnOK7wqZ9G3g7mPGZIFOlT+a3/BrVnl2R9b2OxoSxFdHt+Dj+aq5N/4DlmZ34Kq6/1yGZAKjQlSRMaEvKSaFx3la+jj3P61BMJTAh/gqWRnfkliNv0DRno9fhmACwBGWCpn/GDI5IArOt7T1TDvIkkucT/0qaxPO3w0/Zs6PCgCUoExzpBzgzazbfx5xDpsR5HY2pJA5E1OKFxL/SOHczIw7b9aiKzhKUCY5fPiGWLLsWYMrd0pjOfBx/NX2yZnFxhr8Vh00osgRlAk8VFr/HusgWrIs61etoTCU0Pv5K5kX34Ma0t+D3H70Ox5wgS1Am8LYtgZ2/MMPOnoxHVCJ4IfEvbItsCJ9cDwc2lzyRCTmWoEzg/fweRCfwfcw5XkdiKrH0iASerPoQ5GTBuGsgy27irWgsQZnASt8Py8dD+0vtse7Gc1sjG8Nlb8L2ZfDZzZCX63VIphQsQZnAWvIhZKdBtz96HYkxjtYXQP8n4bcvYObDJY9vQoYlKBM4ebmw4E1o2hManOZ1NMb8T48/QdebYe4o5ztqKgRLUCZwVs+AAxuhu509mRAjAgOeglYDYNp9znfVhDxLUCZwFrwB1RpBm4u8jsSY40VGwWX/gZM6wPjrYdM8ryMyJbAEZQJj12+wfhYk3+D8EBgTimIT4eoJUK0hfPQH2L7c64hMMSxBmcCYOwqi4uD0YV5HYkzxEuvBdZMhtip8eCnsWet1RKYIlqBM2aVug2VjofO1UKWO19EYU7IaTeC6SU6rJ+8Ptht5Q5QlKFN2c18FzYOet3sdiTH+q9MSrp0ImYfg3YGw3x7REWosQZmySdsHi9+F9pdBzWZeR2NM6TQ4zTmTyjjgJKl9v3sdkfFhCcqUzcL/QNZh6HW315EYc2IadYHrP3e+x+8OhL3rvI7IuCxBmROXeRjmvw4t+0H9dl5HY8yJa9ARrv8CcjLgnQth5wqvIzJYgjJlMX80pO2Fs//mdSTGlN1J7WHYl85Nve8MgA0/eR1RpRfUBCUiA0RklYisFZH7CxkeKyLj3OHzRaSZ27+ZiKSLyFL3NTqYcZoTkH4A5rwCrS6AxsleR2NMYNRrCzd+BYn14YNLYIU98NBLQUtQIhIJvApcACQBV4lIUoHRbgT2q+qpwIvA0z7D1qlqJ/d1a7DiNCdo7quQcRD6/N3rSIwJrBpN4YYZTrHf+Ous7T4PBfMMqhuwVlXXq2oWMBYYXGCcwcB7bvcEoK+ISBBjMoFwZC/Mew2ShlijsCY8JdRybuZtfQFM/St8cQ/kZnsdVaUTzATVCPC9+22L26/QcVQ1BzgI1HaHNReRJSLyvYicVdgCROQWEVkkIot2794d2OhN0X541nmkxjkPeB2JMcETkwBXfghn3gWL/uMU+R3Z63VUlUqoVpLYDjRV1c7APcDHIlKt4EiqOkZVk1U1uW7duuUeZKW0ezUsfBO6XA/12ngdjTHBFREJ5z8Gl4yBzQvgzXNgx69eR1VpBDNBbQWa+Hxu7PYrdBwRiQKqA3tVNVNV9wKo6mJgHdAqiLEaf331D4hOgD4Peh2JMeWn45VwwzSnmO+tvs7N6apeRxX2gpmgFgItRaS5iMQAQ4GCVWKmANe73ZcD36qqikhdt5IFInIK0BJYH8RYjT/Wfg1rZkDveyHRzlhNJdPodPjjD9C0B3x+F3x6I2Skeh1VWAtagnKvKd0OzABWAuNVNUVEHhORQe5o/wFqi8hanKK8/KrovYHlIrIUp/LEraq6L1ixGj/kZML0B6Bmc3sgoam8EuvBNRPh3IcgZSK80Ru2LvY6qrAV1Af3qOpUYGqBfg/7dGcAVxQy3afAp8GMzZTSjy/AntVw9acQFet1NMYE1zsDSx6nXjvYswrePBeqNXFaSJcS/vMP/zIw8VUS9mQ5U7Jdv8GPz0OHK6DleV5HY0zAXDxyduED9lzu1/RVSOMmJnFe6gJ+P5jDS1zNehoXPUFRywM+v6OXX8usTCxBmeLl5Trl7bGJ0P/fXkdjTPmo09Kv0Y4AL9OROVkLuP3wKzyvLzIx7lLGJQwlU+KCG2MlEKrVzE2omDMSNs9zkpNVjDCmUAtjujGixuvMiu3DFRmf8NqBW+mROcdq+pWRJShTtG1L4NvHIWkwdBzqdTTGhLTDEVV5OfHP/K3a0xyRKjx4+EkeOfQojXK3eB1ahWUJyhQu6wh8ehNUqQcXveS08GyMKdGK6Pb8ufrLvJVwE0k5Kbx64Db+dHgUNfOsInJpWYIyx1OFKXc4D267ZLTTLpkxxm+5EsXk+Ev4Y403mRZ7IednfsWY/TdxddoHxOeleR1ehWEJyhxvzivw66fQ9yE45WyvozGmwjoQUZM3Em/jTzVGszCmG0PTx/L2gWFcnfYBVfPsJt+SWIIyx1r7NXz9qHPdqdc9XkdjTFjYHtmQZ6rez93VX2Z5dEeGpo/lP/uHc8ORt6iVZw3QFsWqmZv/2boYxl0H9ZJg8Gt23cmYAFsXdSr/rvogTXM2cnn6eAZlTOaijM+ZE3MmbI6Bxl3tuPNhZ1DGsWcNfHQFVKkN13zq3PdkjAmKTVEn80LVe7m1xhimxg0kOXsh/Od8GHMOLP0Ysuw6FViCMuAkp/cGAQLXToKqJ3kdkTGVwo7IBrxV5RaG1XwfBj4P2ekw6TZ4rhVMHgEbK/e9VFbEV9nt+BU+GOIcBNdNhtotvI7ImEonQ+Kh602QfCNs/AmW/hdSJsGSD50Gmk/7A7QdBPXbVaoiQEtQldm67+CT6yG6ipOc6tojt4zxlAg06+W8LnwGVn7uFPn98Cx8/zTUOsVJVEmDoGGXsE9WlqDCjT+tMKvCoe2wf73z8MHaDeGLPxc/jbXCbEz5iqnitODScSgc3g2/fQErp8DcUfDTS5B4ErQ4F07t67yH4f2KlqAqm9xs2LcW0vZCfC2o0woi7GtgTEhLrAvJw51X2j5YPR3WfAWrpsKyjwGBhp3hlHPg5J5ObcD4Gh4HXXb2yxRuijrTUYVV05wzpYxUOO9R6HkXRFg9GWMqlIRa0On/nFdertNm5tpvYN03zk32s18ABOq3d57+27QHNE6GGidXuCJBS1AVUJHPsClCo9wt3HxkDKdnL2Zj5Mm8UPV51i9pAUvmlDkWe4aNMR6KiISZj7jd0dCoK2QegsxUOLgZFqXAwjfd4VFQuyUk1PZv3iFQrG8JKow1ydnEFenj6J31AxkSx5sJN/Nl3EXkiu12YyqSYv+UFvNwxQhyac42WrKJVrkbmbzrHDbS0L+FFrHM8vxTar9UYSZKs+mRNZd+mV/ROXsJ6cQxOW4IE+Mv5UBETa/DM8YEWjEPV8wD1tGGdcD0cgsocCxBhYEqeYfokv0z3bPmc3r2IhL1CLsi6vJR/NVMjRtIakR1r0M0xphSC2qCEpEBwMtAJPCWqj5VYHgs8D5wOrAXuFJVN7jDHgBuBHKBO1V1RjBjrTCO7CUpO4XGuZtpnbOKNjm/0TR3EwAHpDrzYs7gh5izWRbdkTyJ9DhYY4w5cUFLUCISCbwKnA9sARaKyBRVXeEz2o3AflU9VUSGAk8DV4pIEjAUaAc0BL4WkVaqmhuseIMiLw80F/JynNo2muu+5/2vX1620+5W1mH3dQQyDzsXOo/sgsM74fAuOLQDDmyE9P087c4+VaqyKqo138eczfLo01gd1dqSkjEmbATzDKobsFZV1wOIyFhgMOCboAYDj7rdE4BRIiJu/7Gqmgn8LiJr3fnNDWK8xfPnBtiMVNj5CxCgtrMkAqrUhcR6zk15DTtDnVY8OieLLZGN2RlRv8JVGzXGGH8FM0E1Ajb7fN4CdC9qHFXNEZGDQG23/7wC0zYquAARuQW4xf14WERWBSb047WpE9E6GPP9bU9eCTEfANYU7FkH2BOMeEpL7iyXxZR5fcspzoCQO0Nn/xYngNs0qOsbgvu+0PUNwTgLdQJxFra+J/szYYWuJKGqY4AxXsdR3kRkkaomex1HebH1DW+2vuGtLOsbzGYEtgJNfD43dvsVOo6IRAHVcSpL+DOtMcaYMBbMBLUQaCkizUUkBqfSw5QC40wBrne7Lwe+VVV1+w8VkVgRaQ60BBYEMVZjjDEhJmhFfO41pduBGTjVzN9W1RQReQxYpKpTgP8AH7iVIPbhJDHc8cbjVKjIAUZUuBp8wVXZijVtfcObrW94O+H1Fa3ET2s0xhgTuqwpa2OMMSHJEpQxxpiQZAkqRInI2yKyS0R+LWK4iMgrIrJWRJaLSJfyjjGQ/Fjfc0TkoIgsdV8Pl3eMgSQiTUTkOxFZISIpInJXIeOEzT72c33DZh+LSJyILBCRZe76/rOQcWJFZJy7f+eLSLPyjzQw/FzfYSKy22f/3lTijFXVXiH4AnoDXYBfixh+ITANEKAHMN/rmIO8vucAX3gdZwDXtwHQxe2uCqwGksJ1H/u5vmGzj919luh2RwPzgR4FxvkTMNrtHgqM8zruIK/vMGBUaeZrZ1AhSlV/wKnZWJTBwPvqmAfUEJEG5RNd4PmxvmFFVber6s9u9yFgJce3lhI2+9jP9Q0b7j477H6Mdl8Fa6QNBt5zuycAfd2m3iocP9e31CxBVVyFNSUVtge86wy3CGGaiLTzOphAcYt2OuP86/QVlvu4mPWFMNrHIhIpIkuBXcBMVS1y/6pqDpDf1FuF5Mf6AlzmFldPEJEmhQw/hiUoU1H8DJysqh2BkcAkj+MJCBFJBD4F7lbVVK/jCbYS1jes9rGq5qpqJ5yWcLqJSHuvYwomP9b3c6CZqp4GzOR/Z49FsgRVcVWq5qBUNTW/CEFVpwLRIlLH47DKRESicX6sP1LVzwoZJaz2cUnrG477GEBVDwDfAQMKDCqqqbcKraj1VdW96jyhAuAtnOcAFssSVMU1BbjOrenVAzioqtu9DipYROSk/PJ5EemG892tsAezuy7/AVaq6gtFjBY2+9if9Q2nfSwidUWkhtsdj/NcvN8KjFZUU28Vjj/rW+D66SCc65DFqtCtmYczEfkvTq2mOiKyBXgE58IjqjoamIpTy2stkAYM9ybSwPBjfS8HbhORHCAdGFpRD2bXmcC1wC9uuT3A34GmEJb72J/1Dad93AB4T5wHt0YA41X1C/GjqbcKyp/1vVNEBuE0X7cPp1ZfsaypI2OMMSHJiviMMcaEJEtQxhhjQpIlKGOMMSHJEpQxxpiQZAnKGGNMSLIEZYwfRKSxiEwWkTUisk5EXhaRmBKmqSEif/L53FBEJgQ/2mNiGCQi95fnMo0JFKtmbkwJ3JtH5wOvq+o77r0eY4B9qnpvMdM1w2md25MmbkQkym3jzZgKyc6gjCnZuUCGqr4DTptjwJ+BG0QkwX3OzWQRmeWeYT3iTvcU0MJ99s2zItJM3Odduc/PeUdEfhGRJSLSx+0/TEQ+E5Hp7ryeKSwgEdkgIs+40y8QkVPd/u+KyGgRmQ88485vlDusvohMdBtjXSYiPd3+17jzWCoiML9qxAAAAjRJREFUb7gJ2BjPWUsSxpSsHbDYt4eqporIJuBUt1c3oD1Oiw8LReRL4H6gvduAZv4ZVb4Rzmy0g4i0Ab4SkVbusE44rX1nAqtEZKSq+rZqnu+gO/11wEvARW7/xkBPVc0V+f/27t81iiAM4/j3sRRRTGGpYKd1GjGNBnsDSgqbtIKQWvA/EBuL2KoIIqL/gKQQCYI/mlOOqCm01yKIGpHwWMycLIfH5oguqzyfZpdhbmd34Xh5Z4d3tNTofx14bHuhBqF9ko4Bi8BJ2z8krQAXgNvTvKCIvyEZVMSf8agWw/wGPATmWvrPAXcAbK8DH4BRgFq1vWl7CxgCRyZc427jeKLRfr9meeNOAzfqmNu2N4F5StHO57UE0TxwtOXeIzqRDCqi3ZBSJ+4XSfspdeQ2KDsBj3/M3c3H3e+N820m/0894fzLFGMJuGX78hS/iehEMqiIdqvA3jqVRp0euwbctP219jkjaaZWcj4LrAGfKdub/84TylQadWrvMPBmyvtabByf7vA5Lo6eQdKB2nZO0qHaPiNpUsYW0akEqIgWtaL2AnBe0jvgLbBFqcY98oyy19EAeGD7he1PwJqk15Kujl12Bdgj6RVwD1hq7JWzUwclDYBlyqKNNsvAqTrmS+C47SFwhfINbEDZSO6f3FY+/j9ZZh6xS3UhwqztSx2O+b6O+bGrMSO6lgwqIiJ6KRlURET0UjKoiIjopQSoiIjopQSoiIjopQSoiIjopQSoiIjopZ9g7mCAeyMYFgAAAABJRU5ErkJggg==\n",
      "text/plain": [
       "<Figure size 432x288 with 1 Axes>"
      ]
     },
     "metadata": {
      "needs_background": "light"
     },
     "output_type": "display_data"
    }
   ],
   "source": [
    "'''\n",
    "Plot the probability of each unary state on the corresponding option price and\n",
    "check against the precise classical probability drsitribution function for the \n",
    "target log-normal distribution\n",
    "\n",
    "'''\n",
    "width = (S[1]-S[0])/1.2\n",
    "\n",
    "fig, ax = plt.subplots()\n",
    "ax.bar(S, prob_sim, width, label='Quantum', alpha=0.8)\n",
    "ax.scatter(S, ln, s=500, color='C1', label='Exact', marker='_')\n",
    "ax.plot(Sp, lnp* (S[1]- S[0]) / (Sp[1] - Sp[0]), 'C1')\n",
    "plt.ylabel('Probability')\n",
    "plt.xlabel('Option price')\n",
    "plt.title('Option price distribution for {} qubits - qasm_simulator'.format(qu))\n",
    "ax.legend()\n",
    "fig.tight_layout()\n",
    "\n",
    "plt.show()"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The simulation using the calculated exact parameters for the partial-SWAP gate yields a good approximation of the target lognormal distribution. Therefore, we can keep going and apply the expected payoff calculator to this result.\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Expected payoff calculator\n",
    "\n",
    "Once the probability distribution has been succesfully loaded, an ancilla is added to the system\n",
    "$$|\\Psi\\rangle=\\sum_i^n\\sqrt{p_i}|q_i\\rangle|0\\rangle.$$\n",
    "In order to encode the expected payoff into the probability of measuring $|1\\rangle$ in the ancilliary qubit, one must first find the qubits that are mapped to asset prices above the strike $K$, as this example focuses on European _call_ options. In the unary case this is easy, as each qubit is mapped to an asset price. Therefore, the algorithm consists on applying $cR_y$ gates targetting the ancilla, controlled exclusively by qubits with an associated asset price that is higher than the strike $K$. This circuit, for $n$ qubits plus the ancilla, takes the form\n",
    "<img src=\"Images/Payoff_text.png\" width=\"600\" title=\"Payoff circuit\"/>\n",
    "\n",
    "In order to finalize the circuit one must set the angles on the $cR_y$ gates to successfully encode the expected payoff. In this case, each angle will be\n",
    "$$\\theta_i=2\\arcsin{\\sqrt{\\frac{S_i-K}{S_n-K}}},$$\n",
    "where the $S_n-K$ term is set for normalization purposes. After the circuit is applied, the final state of the quantum system will be\n",
    "$$|\\Psi\\rangle=\\sum_{i<K}^n\\sqrt{p_i}|q_i\\rangle|0\\rangle+\\sum_{i>K}^n\\sqrt{p_i}\\cos{\\frac{\\theta_i}{2}}|q_i\\rangle|0\\rangle+\\sum_{i>K}^n\\sqrt{p_i}\\sin{\\frac{\\theta_i}{2}}|q_i\\rangle|1\\rangle.$$\n",
    "Which results in a probability of measuring $|1\\rangle$ in the ancilla of\n",
    "$$P(|1\\rangle) = \\sum_{S_i > K} p_i \\frac{S_i - K}{S_{max} - K},$$\n",
    "where the expected payoff can be extracted by removing the normalization piece."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 20,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Circuit that codifies the expected payoff of the option into the probability of a\n",
    "single ancilla qubit\n",
    "\n",
    "'''\n",
    "def payoff_circuit(qubits, K, S):\n",
    "    C = QuantumCircuit(qubits+1)\n",
    "    for i in range(qubits): #Determine the first qubit's price that \n",
    "        qK = i              #surpasses the strike price\n",
    "        if K<S[i]:\n",
    "            break\n",
    "    C.barrier(range(qubits+1))\n",
    "    for i in range(qK, qubits):                              #Control-RY rotations controled by states\n",
    "        angle = 2 * np.arcsin(np.sqrt((S[i]-K)/(S[qubits-1]-K))) #with higher value than the strike\n",
    "        C.cu3(angle, 0, 0, i, qubits)                        #targeting the ancilla qubit\n",
    "    return C"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Error mitigation\n",
    "\n",
    "Here is where the post-selection scheme appears. As the amplitude distributor works solely in a unary representation, if the read-out of any measurement of the first $n$ qubits contains more than one $|1\\rangle$, there has been an error in the process. Therefore, that run can be discarded.\n",
    "\n",
    "This simple post-selection scheme does, in fact, reduce the error of the outcome. \n"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 21,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Circuit that measures the ancilla qubit to sample the expected payoff of the option\n",
    "\n",
    "'''\n",
    "def measure_payoff(qubits):\n",
    "    '''\n",
    "    Measuements needed for the post-selection step of the algorithm\n",
    "    \n",
    "    '''\n",
    "    C = QuantumCircuit(qubits+1, qubits+1)\n",
    "    C.measure(range(qubits+1), range(qubits+1))\n",
    "    return C\n"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Build the complete circuit\n",
    "\n",
    "Now the complete circuit for expected payoff calculation can be concatenated together. The final circuit includes the probability loading algorithm for the unary representation and the payoff circuit. An example for a complete circuit for 7 qubits plus an ancilla.\n",
    "<img src=\"Images/Total_text.png\" width=\"800\" title=\"Total circuit\"/>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 22,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Total circuit containing the probability loading and the payoff estimation measuring\n",
    "the ancilla qubit.\n",
    "\n",
    "'''\n",
    "circ_payoff_sim = prob_loading + payoff_circuit(qu, K, S) + measure_payoff(qu)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "###### Simulation"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 23,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Run the complete payoff expectation circuit through a simulator \n",
    "and sample the ancilla qubit\n",
    "\n",
    "'''\n",
    "job_payoff_sim = execute(circ_payoff_sim, backend_sim, shots=samples)\n",
    "result_payoff_sim = job_payoff_sim.result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 24,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_payoff_sim = result_payoff_sim.get_counts(circ_payoff_sim)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 25,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Post-selection step where states that do not belong to the unary representation of\n",
    "the asset value are discarded.\n",
    "\n",
    "'''\n",
    "ones=0\n",
    "zeroes=0\n",
    "for key in counts_payoff_sim.keys(): #Post-selection step\n",
    "    unary = 0\n",
    "    for i in range(1,qu+1):\n",
    "        unary+=int(key[i])\n",
    "    if unary==1:\n",
    "        if int(key[0])==0:\n",
    "            zeroes+=counts_payoff_sim.get(key)\n",
    "        else:\n",
    "            ones+=counts_payoff_sim.get(key)\n",
    "\n",
    "qu_payoff_sim = ones * (S[qu - 1]-K) / (ones+zeroes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 26,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Calculation of the expected payoff using the precise probability distribution\n",
    "computed classically\n",
    "\n",
    "'''\n",
    "cl_payoff = 0\n",
    "for i in range(len(Sp)):\n",
    "    if K<Sp[i]:\n",
    "        cl_payoff += lnp[i]*(Sp[i]-K)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 27,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With precision 10000\n",
      "Classical Payoff: 0.21755186788182482\n",
      "\n",
      "With 8 qubits and 10000 samples\n",
      "Quantum Payoff: 0.22465089021855428\n",
      "\n",
      "Percentage off: 3.2631401448530344%\n"
     ]
    }
   ],
   "source": [
    "print('With precision {}'.format(samples))\n",
    "print('Classical Payoff: {}'.format(cl_payoff))\n",
    "print('')\n",
    "print('With {} qubits and {} samples'.format(qu, samples))\n",
    "print('Quantum Payoff: {}'.format(qu_payoff_sim))\n",
    "print('')\n",
    "print('Percentage off: {}%'.format(np.abs(100*(cl_payoff-qu_payoff_sim)/cl_payoff)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "___Results:___ \n",
    "\n",
    "Even with an exponential loss of precision, it seems that the final result for the expected payoff does not differ much from the precise value computed by classical means.\n",
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "In order to test the resiliance of the algorithm under noise, we could run it through a real quantum device. But if we want more control over the errors, we can define a noise model in order to run noisy simulations. "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 28,
   "metadata": {},
   "outputs": [],
   "source": [
    "from qiskit import QuantumCircuit, execute, Aer\n",
    "from qiskit.providers.aer.noise.errors import pauli_error\n",
    "from qiskit.providers.aer.noise import NoiseModel"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 29,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Build a noise model with tunable bitflip, phaseflip and measurement errors.\n",
    "\n",
    "'''\n",
    "def noise_model_total(qubits, t_bit, t_phase, t_meas):\n",
    "    noise_model = NoiseModel()\n",
    "    bit_error = t_bit\n",
    "    phase_error = t_phase\n",
    "    cnotb_error = 2*t_bit\n",
    "    cnotp_error = 2*t_phase\n",
    "    p_measure = t_meas\n",
    "    bit_flip = pauli_error([('X', bit_error), ('I', 1 - bit_error)])\n",
    "    phase_flip = pauli_error([('Z', phase_error), ('I', 1 - phase_error)])\n",
    "    bitphase_flip = bit_flip.compose(phase_flip)\n",
    "    cnot_bit = pauli_error([('X', cnotb_error), ('I', 1 - cnotb_error)])\n",
    "    cnot_phase = pauli_error([('Z', cnotp_error), ('I', 1 - cnotp_error)])\n",
    "    cnot_error = cnot_bit.compose(cnot_phase)\n",
    "    cnot_error = cnot_error.tensor(cnot_error)\n",
    "    noise_model.add_all_qubit_quantum_error(cnot_error, ['cx'], warnings=False)\n",
    "    noise_model.add_all_qubit_quantum_error(bitphase_flip, [\"u1\", \"u2\", \"u3\"])\n",
    "    measure_error = pauli_error([('X',p_measure), ('I', 1 - p_measure)])\n",
    "    noise_model.add_all_qubit_quantum_error(measure_error, \"measure\")\n",
    "\n",
    "    return noise_model"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Noisy simulation "
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 30,
   "metadata": {},
   "outputs": [],
   "source": [
    "circ_payoff_noise = prob_loading + payoff_circuit(qu, K, S) + measure_payoff(qu)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 31,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Determine the noise model we want to use.\n",
    "\n",
    "In this example: 0.3% bitflip and phaseflip errors\n",
    "                 1% measurement error\n",
    "\n",
    "'''\n",
    "t_bit = 0.003\n",
    "t_phase = 0.003\n",
    "t_meas = 0.01\n",
    "\n",
    "noise_model = noise_model_total(qu, t_bit, t_phase, t_meas)\n",
    "basis_gates = noise_model.basis_gates"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 32,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Run the complete payoff expectation circuit through noisy simulation\n",
    "following the determined noise model\n",
    "\n",
    "'''\n",
    "backend_noise = Aer.get_backend('qasm_simulator')\n",
    "job_payoff_noise = execute(circ_payoff_noise, backend_noise, shots=samples, basis_gates=basis_gates, noise_model=noise_model)\n",
    "result_payoff_noise = job_payoff_noise.result()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 33,
   "metadata": {},
   "outputs": [],
   "source": [
    "counts_payoff_noise = result_payoff_noise.get_counts(circ_payoff_noise)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 34,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Post-selection step where states that do not belong to the unary representation of\n",
    "the asset value are discarded.\n",
    "\n",
    "'''\n",
    "ones=0\n",
    "zeroes=0\n",
    "\n",
    "for key in counts_payoff_noise.keys():\n",
    "    unary = 0\n",
    "    for i in range(1,qu+1):\n",
    "        unary+=int(key[i])\n",
    "    if unary==1:\n",
    "        if int(key[0])==0:\n",
    "            zeroes+=counts_payoff_noise.get(key)\n",
    "        else:\n",
    "            ones+=counts_payoff_noise.get(key)\n",
    "\n",
    "qu_payoff_post = ones * (S[qu - 1]-K) / (ones+zeroes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 35,
   "metadata": {},
   "outputs": [],
   "source": [
    "'''\n",
    "Result without the post-selection scheme\n",
    "\n",
    "'''\n",
    "ones=0\n",
    "zeroes=0\n",
    "\n",
    "for key in counts_payoff_noise.keys():\n",
    "    if int(key[0])==0:\n",
    "        zeroes+=counts_payoff_noise.get(key)\n",
    "    else:\n",
    "        ones+=counts_payoff_noise.get(key)\n",
    "\n",
    "qu_payoff_no_post = ones * (S[qu - 1]-K) / (ones+zeroes)"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 36,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "With precision 10000\n",
      "Classical Pay-off: 0.21755186788182482\n",
      "\n",
      "Noise simulation with 0.003% bitflip, 0.003% phaseflip and 0.01% measurement errors:\n",
      "With precision 8 qubits and 10000 samples\n",
      "\n",
      "With post-processing\n",
      "\n",
      "Quantum Pay-off: 0.2937883678863928\n",
      "\n",
      "Percentage off: 35.04290758191972%\n",
      "\n",
      "Without post-processing\n",
      "\n",
      "Quantum Pay-off: 0.3646350391506456\n",
      "\n",
      "Percentage off: 67.60832379923166%\n"
     ]
    }
   ],
   "source": [
    "print('With precision {}'.format(samples))\n",
    "print('Classical Pay-off: {}'.format(cl_payoff))\n",
    "print('')\n",
    "print('Noise simulation with {}% bitflip, {}% phaseflip and {}% measurement errors:'.format(t_bit, t_phase, t_meas))\n",
    "print('With precision {} qubits and {} samples'.format(qu, samples))\n",
    "print('')\n",
    "print('With post-processing')\n",
    "print('')\n",
    "print('Quantum Pay-off: {}'.format(qu_payoff_post))\n",
    "print('')\n",
    "print('Percentage off: {}%'.format(np.abs(100*(cl_payoff-qu_payoff_post)/cl_payoff)))\n",
    "print('')\n",
    "print('Without post-processing')\n",
    "print('')\n",
    "print('Quantum Pay-off: {}'.format(qu_payoff_no_post))\n",
    "print('')\n",
    "print('Percentage off: {}%'.format(np.abs(100*(cl_payoff-qu_payoff_no_post)/cl_payoff)))"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "The post-processing scheme reduced over 30% of the error found in this highly noisy simulation."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "***"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "#### Summary\n",
    "\n",
    "The unary representation may need a larger amount of qubits to reach the same precision than the binary scheme allows for. However, the simplicity of the unary circuit and the ability to post-select the outcomes to achieve error mitigation could make it ideal for NISQ devices."
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### _References_\n",
    "\n",
    "[1] R. Orus, S. Mugel and E. Lizaso, *Quantum computing for finance: overview and prospects*, Reviews in Physics 4 (2019) 100028, [arXiv:1807.03890](https://arxiv.org/abs/1807.03890). <br>\n",
    "\n",
    "[2] P. Rebentrost, B. Gupt and T. R. Bromley, *Quantum computational finance: Monte Carlo pricing of financial derivatives*, Phys. Rev. A 98, 022321 (2018), [arXiv:1805.00109](https://arxiv.org/abs/1805.00109). <br>\n",
    "\n",
    "[3] N. Stamatopoulos, D. J. Egger, Y. Sun, C. Zoufal, R. Iten, N. Shen and S. Woerner, *Option Pricing using Quantum Computers* (2019), [arXiv:1905.02666](https://arxiv.org/abs/1905.02666). <br>\n",
    "\n",
    "[4] F. Black and M. Scholes, *The Pricing of Options and Corporate Liabilities*,  The Journal of Political Economy, Vol. 81, No. 3, pp. 637-654 (1973), [JSTOR](http://www.jstor.org/stable/1831029). <br>\n",
    "\n",
    "[5] L. Grover and T. Rudolph, *Creating superpositions that correspond to efficiently integrable probability distributions* (2002), [arXiv:quant-ph/0208112](https://arxiv.org/abs/quant-ph/0208112). <br>\n",
    "\n",
    "[6] S. A. Cuccaro, T. G. Draper, S. A. Kutin and D. Petrie Moulton, *A new quantum ripple-carry addition circuit* (2004), [arXiv:quant-ph/0410184](https://arxiv.org/abs/quant-ph/0410184). <br>\n",
    "\n",
    "[7] M.A. Nielsen and I.L. Chuang, *Quantum computation and quantum information*, Cambridge UniversityPress (2010). <br>"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 37,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/html": [
       "<h3>Version Information</h3><table><tr><th>Qiskit Software</th><th>Version</th></tr><tr><td>Qiskit</td><td>0.13.0</td></tr><tr><td>Terra</td><td>0.10.0</td></tr><tr><td>Aer</td><td>0.3.2</td></tr><tr><td>Ignis</td><td>0.2.0</td></tr><tr><td>Aqua</td><td>0.6.1</td></tr><tr><td>IBM Q Provider</td><td>0.3.3</td></tr><tr><th>System information</th></tr><tr><td>Python</td><td>3.6.8 (default, Aug 20 2019, 17:12:48) \n",
       "[GCC 8.3.0]</td></tr><tr><td>OS</td><td>Linux</td></tr><tr><td>CPUs</td><td>4</td></tr><tr><td>Memory (Gb)</td><td>31.35873794555664</td></tr><tr><td colspan='2'>Thu Dec 12 17:28:05 2019 CET</td></tr></table>"
      ],
      "text/plain": [
       "<IPython.core.display.HTML object>"
      ]
     },
     "metadata": {},
     "output_type": "display_data"
    }
   ],
   "source": [
    "import qiskit.tools.jupyter\n",
    "%qiskit_version_table"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": null,
   "metadata": {},
   "outputs": [],
   "source": []
  }
 ],
 "metadata": {
  "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.0"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}