diff --git a/tnbs/BTC_02_AE/QQuantLib/NEASQC_Benchmark_AE_usingintegrals.ipynb b/tnbs/BTC_02_AE/QQuantLib/notebooks/01_AmplitudeEstimationKernel.ipynb similarity index 57% rename from tnbs/BTC_02_AE/QQuantLib/NEASQC_Benchmark_AE_usingintegrals.ipynb rename to tnbs/BTC_02_AE/QQuantLib/notebooks/01_AmplitudeEstimationKernel.ipynb index 0813d3d..2351818 100644 --- a/tnbs/BTC_02_AE/QQuantLib/NEASQC_Benchmark_AE_usingintegrals.ipynb +++ b/tnbs/BTC_02_AE/QQuantLib/notebooks/01_AmplitudeEstimationKernel.ipynb @@ -5,19 +5,9 @@ "id": "0618c8ae", "metadata": {}, "source": [ - "# Integral Computation using Amplitude Estimation" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "7071be0e", - "metadata": {}, - "outputs": [], - "source": [ - "import numpy as np\n", - "import matplotlib.pyplot as plt\n", - "%matplotlib inline" + "# The Amplitude Estimation Kernel\n", + "\n", + "In this notebook the **Amplitude Estimation Kernel** is introduced. " ] }, { @@ -134,1222 +124,34 @@ "## 5. Amplitude Estimation without Phase Estimation\n", "\n", "\n", - "**Phase Estimation** algorithm is computationally expensive and present quantum computers do not have enough quality, neither qubits, for properly implemented it. There exists, however, several algorithms that can solve the *Amplitude Estimation* problem without the use of **Phase Estimation** and with scaling errors between **MonteCarlo** and **Phase estimation**, this is \n", + "**Phase Estimation** algorithm is computationally expensive and present quantum computers do not have enough quality, nor qubits, to properly implement it. There exists, however, several algorithms that can solve the *Amplitude Estimation* problem without the use of **Phase Estimation** and with scaling errors between **Monte-Carlo** and **Phase estimation**, this is \n", "\n", "$$\\frac{1}{N} < \\epsilon < \\frac{1}{\\sqrt{N}} $$\n", "\n", - "These algorithms are known as *Amplitude Estimation* algorithms. Main idea of these algorithms is taking advantage of the fact:\n", + "These algorithms are known as *Amplitude Estimation* algorithms. The main idea of these algorithms is to take advantage of the facts:\n", "\n", "$$\\mathbf{G}^k |\\Psi\\rangle = \\mathbf{G}^k \\mathbf{A} |0\\rangle_n= \\sin\\big((2k+1)\\theta\\big) |\\Psi_0\\rangle + \\cos\\big((2k+1)\\theta\\big)|\\Psi_1\\rangle$$\n", "\n", "and in the use of very smart strategies for selecting $k$ in order to maximize the probability of measuring the $|\\Psi_0\\rangle$:\n", "\n", - "$$P[|\\Psi_0\\rangle] = \\sin^2\\big((2k+1)\\theta\\big)$$" - ] - }, - { - "cell_type": "markdown", - "id": "c94e6e4c", - "metadata": {}, - "source": [ - "## 6. Benchmarking Amplitude Estimation algorithms.\n", - "\n", - "Here we propose a benchmark procedure for evaluating *Amplitude Estimation* algorithms by computing integrals. " - ] - }, - { - "cell_type": "markdown", - "id": "011635f6", - "metadata": {}, - "source": [ - "### 6.1. Description of the problem\n", - "\n", - "The benchmark problem is the computation of the integral of a function $f(x)$ in a closed interval $[a,b] \\subset \\mathbf{R}$:\n", - "\n", - "$$I = \\int_a^bf(x)dx$$\n", - "\n", - "In particular we propose to use $f(x) = \\sin x$, whose integral is very well known:\n", - "\n", - "$$I = \\int_a^{b}\\sin(x)dx = -\\cos x |_a^b = \\cos(a)-\\cos(b)$$\n", - "\n", - "Additionally 3 different integration intervals will be used:\n", - "\n", - "1. $[0, \\frac{3\\pi}{8}]$ : $I = \\int_0^{\\frac{3\\pi}{8}}\\sin(x)dx = 0.6173165676349102$\n", - "2. $[\\frac{3\\pi}{4}, \\frac{9\\pi}{8}]$ : $I = \\int_{\\frac{3\\pi}{4}}^{\\frac{9\\pi}{8}}\\sin(x)dx = 0.2928932188134523$\n", - "3. $[\\pi, \\frac{5\\pi}{4}]$ : $I = \\int_{\\pi}^{\\frac{5\\pi}{4}}\\sin(x)dx = -0.2928932188134523$\n", - "\n", - "So the *Amplitude Estimation* algorithm should be used for computing each on of the 3 integrals. \n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "871a0510", - "metadata": {}, - "outputs": [], - "source": [ - "def sin_integral(a,b):\n", - " return np.cos(a)-np.cos(b)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c7aec8ad", - "metadata": {}, - "outputs": [], - "source": [ - "start = [0.0, 3.0*np.pi/4.0, np.pi]\n", - "end = [3.0*np.pi/8.0, 9.0*np.pi/8.0, 5.0*np.pi/4.0]" - ] - }, - { - "cell_type": "markdown", - "id": "148bbae2", - "metadata": {}, - "source": [ - "#### $1^{st}$ Domain: $[0, \\frac{3\\pi}{8}]$ \n", - "\n", - "In this domain the function is strictly posiitve defined." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "29114c8f", - "metadata": { - "scrolled": false - }, - "outputs": [], - "source": [ - "a = start[0]\n", - "b = end[0]\n", - "domain = np.linspace(a,b, 100)\n", - "plt.plot(domain, np.sin(domain))\n", - "#plt.show()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "d8fcfba3", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "print('Integral in first domain: {}'.format(sin_integral(a,b)))" - ] - }, - { - "cell_type": "markdown", - "id": "0670f881", - "metadata": {}, - "source": [ - "#### $2^{nd}$ Domain: $[\\frac{3\\pi}{4}, \\frac{9\\pi}{8}]$ \n", - "\n", - "In this domain the function is not strictly positive defined but its integral is positive." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "47ed9e21", - "metadata": {}, - "outputs": [], - "source": [ - "a = start[1]\n", - "b = end[1]\n", - "domain = np.linspace(a,b, 100)\n", - "plt.plot(domain, np.sin(domain))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cf914ead", - "metadata": {}, - "outputs": [], - "source": [ - "print('Integral in second domain: {}'.format(sin_integral(a,b)))" - ] - }, - { - "cell_type": "markdown", - "id": "387e2fd6", - "metadata": {}, - "source": [ - "#### $3^{rd}$ Domain: $[\\pi, \\frac{5\\pi}{4}]$\n", - "\n", - "In this domain the function is not strictly positve defined and the intgral is negative!" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "9bc7ab74", - "metadata": {}, - "outputs": [], - "source": [ - "a = start[1]\n", - "b = end[1]\n", - "domain = np.linspace(a,b, 100)\n", - "plt.plot(domain, np.sin(domain))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b544a0bc", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "print('Integral in thrid domain: {}'.format(sin_integral(a,b)))" - ] - }, - { - "cell_type": "markdown", - "id": "4725ec35", - "metadata": {}, - "source": [ - "### 6.2. Domain Discretization\n", - "\n", - "First thing to do for computing the integral in a computer is the discretization of the domain. In the benchmark we allways discretize the domain in $2^n$ intervals, with $n \\in \\mathbf{N}$:\n", - "\n", - "$$\\{[x_0, x_1], [x_1, x_2], ..., [x_{2^n-1}, x_{2^n}]\\}$$ \n", - "\n", - "Where\n", - "\n", - "1. $x_{i+1} < x_{i}$\n", - "2. $a = x_0$\n", - "3. $b = x_{2^n}$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c4d9ce5f", - "metadata": {}, - "outputs": [], - "source": [ - "#First integration domain\n", - "a = start[0]\n", - "b = end[0]\n", - "#For fix the number of discretization intervals\n", - "n=4\n", - "domain_x = np.linspace(a, b, 2**n+1)" - ] - }, - { - "cell_type": "markdown", - "id": "8f43650e", - "metadata": {}, - "source": [ - "### 6.3. Function discretization\n", - "\n", - "Now using the domain discretization we need to discretized $f(x)$. In our procedure following arrays should be constructed:\n", - "\n", - "1. $\\Delta x_i = x_{i+1} - x_{i} = \\frac{b-a}{2^n}$\n", - "2. $f_{x_i} = \\frac{f(x_{i+1}) + f(x_{i})}{2}$\n", - "3. $f_{x_i} \\Delta x_i = f_{x_i} \\frac{b-a}{2^n}$\n", - "\n", - "Using these computed arrays the desired integral can be approximated by Riemann sum:\n", - "\n", - "$$S_{[a,b]} = \\sum_{i=0}^{2^n-1} f_{x_i} \\Delta x_i$$\n", - "\n", - "When $\\Delta x_i \\rightarrow 0$ then $I =\\int_a^{b}\\sin(x)dx \\approx S_{[a,b]}$.\n", - "\n", - "Using $\\Delta x_i = \\frac{b-a}{2^n}$ then we can write down:\n", - "\n", - "\n", - "$$S_{[a,b]} = \\sum_{i=0}^{2^n-1} f_{x_i} \\frac{b-a}{2^n} = \\frac{b-a}{2^n} \\sum_{i=0}^{2^n-1} f_{x_i}$$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a30668da", - "metadata": {}, - "outputs": [], - "source": [ - "#The selected fucntion\n", - "f = np.sin\n", - "\n", - "#Discretization of the selected function\n", - "f_x = []\n", - "x_ = []\n", - "for i in range(1, len(domain_x)):\n", - " step_f = (f(domain_x[i]) + f(domain_x[i-1]))/2.0\n", - " #print(i)\n", - " f_x.append(step_f)\n", - " x_.append((domain_x[i] + domain_x[i-1])/2.0)\n", - "f_x = np.array(f_x)\n", - "x_ = np.array(x_)\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "39cabf6d", - "metadata": {}, - "outputs": [], - "source": [ - "plt.plot(domain_x, f(domain_x), '-o')\n", - "plt.plot(x_, f_x, 'o')\n", - "plt.legend(['original', 'half'])" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "0c6a644b", - "metadata": {}, - "outputs": [], - "source": [ - "Riemann = (np.sum(f_x)*(b-a))/2**n\n", - "print(\"Riemann sum integral: {}\".format(Riemann))\n", - "print('Exact Integral in first domain: {}'.format(sin_integral(a,b)))" - ] - }, - { - "cell_type": "markdown", - "id": "785e8e08", - "metadata": {}, - "source": [ - "### 6.4 Array Normalisation\n", - "\n", - "The idea is encode the $2^n$ discretized array $f_{x_i}$ in a $n+1$ qubit circuit. Before doing that we are going to normalise the arrray in the following way:\n", - "\n", - "\n", - "$$f\\_norm_{x_i} = \\frac{f_{x_i}}{\\max(|f_{x_i}|)}$$\n", - "\n", - "If $\\max{|f_{x_i}|} \\leq 1$ then this step can be omitted." - ] - }, - { - "cell_type": "markdown", - "id": "e6f6468f", - "metadata": {}, - "source": [ - "$$S_{[a,b]} = \\frac{b-a}{2^n} \\sum_{i=0}^{2^n-1} f_{x_i} = \\frac{b-a}{2^n} \\sum_{i=0}^{2^n-1} \\max(|f_{x_i}|) f\\_norm_{x_i} =\\frac{\\max(|f_{x_i}|)(b-a)}{2^n} \\sum_{i=0}^{2^n-1} f\\_norm_{x_i}$$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "092a27f7", - "metadata": {}, - "outputs": [], - "source": [ - "normalization = np.max(np.abs(f_x))\n", - "print(\"Normalization constant: {}\".format(normalization))\n", - "#normalization = 1.0\n", - "f_norm_x = f_x/normalization\n", - "Riemann = normalization * np.sum(f_norm_x)*(b-a)/2**n\n", - "#Now we need to be aware of the normalization constant when computing Rieman sum\n", - "print(\"Riemman sum integral: {}\".format(Riemann))\n", - "print('Exact Integral in first domain: {}'.format(sin_integral(a,b)))" - ] - }, - { - "cell_type": "markdown", - "id": "8dc400f5", - "metadata": {}, - "source": [ - "### 6.5 Encoding function in a quantum circuit.\n", - "\n", - "Now we need to codified the $f\\_norm_{x_i}$ array in a quantum circuit.\n", - "\n", - "Several procedures can be used. In this benchmark we are going to use following one:\n", - "\n", - "1. Initialize $n+1$ qubits, where $n$ should be equal to the $n$ of the $2^n$ discretization intervals:\n", - "\n", - "$$|0\\rangle\\otimes|0\\rangle_n \\tag{1}$$\n", - "\n", - "2. Apply the uniform distribution over the first $n$ qubits:\n", - "\n", - "$$\\big(I \\otimes H^{\\otimes n}\\big)\\big(|0\\rangle \\otimes|0\\rangle_{n}\\big) = |0\\rangle \\otimes H^{\\otimes n}|0\\rangle_{n}=\n", - "\\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1}|0\\rangle \\otimes|i\\rangle_{n} \\tag{2}$$\n", - "\n", - "3. Creates an operator $\\mathbf{U}_f$ for encoding the $f\\_norm_{x_i}$. This operator acts in the following way:\n", - "\n", - "$$\\mathbf{U}_f|i\\rangle_n \\otimes |0\\rangle = |i\\rangle_n \\otimes \\big(f\\_norm_{x_i}|0\\rangle + \\beta_i |1\\rangle \\big) \\tag{3}$$\n", - "\n", - "4. In the equation $(3)$ the coeficient of $|1\\rangle$ it is not important for us the only important coeficient it is the $|0\\rangle$ one. \n", - "\n", - "5. Apply the $\\mathbf{U}_f$ operator over $n+1$ qubits:\n", - "\n", - "$$\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n}\\right)|0\\rangle\\otimes|0\\rangle_{n} \\tag{4}$$\n", - "\n", - "6. Applying equation $(2)$ on $(4)$:\n", - "\n", - "$$\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\mathbf{U}_f \\left(\\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} |0\\rangle\\otimes|i\\rangle_{n}\\right) = \\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} \\mathbf{U}_f \\left(|0\\rangle\\otimes|i\\rangle_{n}\\right) \\tag{5}$$\n", - "\n", - "7. Applying equation $(3)$ on $(5)$:\n", - "\n", - "$$\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} |i\\rangle_{n} \\otimes \\left(f\\_norm_{x_i}|0\\rangle + \\beta_i|1\\rangle \\right) \\tag{6}$$\n", - "\n", - "8. Finally the uniform distribution will be applied over the first n qubits again:\n", - "\n", - "$$|\\Psi \\rangle = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} \\tag{7}$$\n", - "\n", - "9. So applying $6$ on $(7)$:\n", - "\n", - "$$|\\Psi \\rangle = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} H^{\\otimes n}|i\\rangle_{n} \\otimes \\left(f\\_norm_{x_i}|0\\rangle + \\beta_i|1\\rangle \\right) \\tag{6}$$\n", - "\n", - "10. We are interested only in $|0\\rangle \\otimes |i\\rangle_{n}$ so we don't need to take into acount other terms, so $(6)$ can be expresed as:\n", - "\n", - "$$|\\Psi \\rangle = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} f\\_norm_{x_i} H^{\\otimes n}|i\\rangle_{n} \\otimes |0\\rangle + \\cdots \\tag{7}$$\n", - "\n", - "11. We know that:\n", - "\n", - "$$H^{\\otimes n} = \\frac{1}{\\sqrt{2^n}} \\sum_{j=0}^{2^n}\\sum_{k=0}^{2^n} (-1)^{jk} |j\\rangle_n {}_{n} \\langle k|$$ \n", - "\n", - "12. So:\n", - "$$H^{\\otimes n} |i\\rangle_n = \\frac{1}{\\sqrt{2^n}} \\sum_{j=0}^{2^n}\\sum_{k=0}^{2^n} (-1)^{jk} |j\\rangle_n {}_{n} \\langle k|i\\rangle_n = \\frac{1}{\\sqrt{2^n}} \\sum_{j=0}^{2^n} (-1)^{ji} |j\\rangle = \\frac{1}{\\sqrt{2^n}} |0\\rangle_n + \\frac{1}{\\sqrt{2^n}} \\sum_{j=1}^{2^n} (-1)^{ji} |j\\rangle \\tag{8}$$ \n", - "\n", - "13. Finally applying $(8)$ in $(7)$ and taking only into account the $|0\\rangle \\otimes |0\\rangle_{n}$ \n", - "\n", - "$$|\\Psi \\rangle = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\frac{1}{2^n} \\sum_{i=0}^{2^{n}-1} f\\_norm_{x_i} |0\\rangle\\otimes|0\\rangle_{n} + \\cdots \\tag{9}$$\n", - "\n", - "14. Now we measured the probability of obtain the state $|\\Psi_0\\rangle = |0\\rangle \\otimes |0\\rangle_n$\n", - "\n", - "$$\\mathbf{P}[|\\Psi_0\\rangle] = \\left| \\; \\langle \\Psi_0\\ |\\Psi\\rangle \\; \\right|^2 = \\left| \\; \\langle \\Psi_0\\ | \\frac{1}{2^n} \\sum_{i=0}^{2^{n}-1} f\\_norm_{x_i} |\\Psi_0\\rangle\\; \\right|^2 = \\left| \\frac{1}{2^n} \\sum_{i=0}^{2^{n}-1} f\\_norm_{x_i} \\right|^2 \\tag{10}$$ \n", - "\n", - "\n", - "15. So reorganising $(10)$:\n", - "\n", - "$$ \\left| \\sum_{i=0}^{2^{n}-1} f\\_norm_{x_i} \\right| = 2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]} \\tag{11}$$\n", - "\n", - "16. Undoing the array normalisation:\n", - "\n", - "$$ \\left| \\sum_{i=0}^{2^{n}-1} \\frac{f_{x_i}}{\\max(f_{x_i})} \\right| = 2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]} \\tag{12}$$\n", - "\n", - "Finally we reorganise $(12)$ and the desired Riemann sum can be obtained:\n", - "\n", - "\n", - "\n", - "$$\\left| \\sum_{i=0}^{2^{n}-1} f_{x_i} \\right| = \\max(f_{x_i}) 2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}$$\n", - "\n", - "As explained before the Riemann sum is:\n", - "\n", - "$$S_{[a,b]} = \\frac{b-a}{2^n} \\sum_{i=0}^{2^n-1} f_{x_i}$$\n", - "\n", - "Then:\n", - "\n", - "$$\\frac{2^n}{b-a} S_{[a,b]} = \\max(f_{x_i}) 2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}$$\n", - "\n", - "So we have obtained an estimation of the desired integral:\n", - "\n", - "$$\\tilde{S}_{[a,b]} = \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}\\right)$$\n", - "\n", - "\n", - "The $2^n$ terms can be erased from equation by for now we are going to keep it!!!\n", - "\n", - "So in our procedure the operator $\\mathbf{A}$ will be:\n", - "\n", - "$$\\mathbf{A}(g_{x_i}) = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)$$\n", - "\n", - "So for each of the 3 integrals an operator $\\mathbf{A}$ should be built.\n", + "$$P[|\\Psi_0\\rangle] = \\sin^2\\big((2k+1)\\theta\\big)$$\n", "\n" ] }, { "cell_type": "markdown", - "id": "bf0e3a5d", - "metadata": {}, - "source": [ - "#### Operator $\\mathbf{U}_f$\n", - "\n", - "Here we are going to explain how to build the operator $\\mathbf{U}_f$:\n", - "\n", - "1. Given the array $f\\_norm_{x_i}$ we need to compute: $\\phi_{x_i} = \\arccos({f\\_norm_{x_i}})$.\n", - "2. for a given state $|i\\rangle_n \\otimes |0\\rangle$ we need to implement a rotation around the *y-axis* over the last qubit ($|0\\rangle$) controlled by the state $|i\\rangle_n$ of $2*\\phi_{x_i}$. So we need to build following operation:\n", - "\n", - "$$|i\\rangle_n \\otimes |0\\rangle \\rightarrow |i\\rangle_n \\otimes \\mathbf{R}_y(2*\\phi_{x_i})|0\\rangle = |i\\rangle_n \\otimes \\left( \\cos(\\phi_{x_i})|0\\rangle + \\sin(\\phi_{x_i})|1\\rangle \\right) $$\n", - "\n", - "3. Now undoing the $\\phi_{x_i}$ and doing $\\beta_i = \\sin(\\phi_{x_i})$ we can obtain the desired operator $\\mathbf{U}_f$:\n", - "\n", - "$$|i\\rangle_n \\otimes \\left(f\\_norm_{x_i} |0\\rangle + \\beta_i|1\\rangle \\right) = \\mathbf{U}_f |i\\rangle_n \\otimes |0\\rangle$$" - ] - }, - { - "cell_type": "markdown", - "id": "ca85032e", - "metadata": {}, - "source": [ - "In the case of **QQuantLib** the complete encoding procedure can be executed in a transparent way by using the python class *Encoding* from *QQuantLib.DL.encoding_protocols* module.\n", - "\n", - "This class creates the operator $\\mathbf{A}$ under the *oracle* property of the *Encoding* class" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a01a0395", - "metadata": {}, - "outputs": [], - "source": [ - "import sys\n", - "sys.path.append(\"../\")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "193f009a", - "metadata": {}, - "outputs": [], - "source": [ - "from QQuantLib.DL.encoding_protocols import Encoding" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "49386077", - "metadata": {}, - "outputs": [], - "source": [ - "encoding_object = Encoding(\n", - " array_function=f_norm_x, \n", - " array_probability=None, \n", - " encoding=2\n", - ")\n", - "encoding_object.run()" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "03299414", - "metadata": {}, - "outputs": [], - "source": [ - "operator_A = encoding_object.oracle\n", - "%qatdisplay operator_A --svg" - ] - }, - { - "cell_type": "markdown", - "id": "440db60f", - "metadata": {}, - "source": [ - "Additionally the *function_gate* attribute give us acces to the operator: $\\mathbf{U}_g$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a08cc3b0", - "metadata": {}, - "outputs": [], - "source": [ - "operator_Uf = encoding_object.function_gate\n", - "%qatdisplay operator_Uf --depth --svg" - ] - }, - { - "cell_type": "markdown", - "id": "36788404", - "metadata": {}, - "source": [ - "## 7. Grover-like operator of $\\mathbf{A}$\n", - "\n", - "For each one of the 3 $\\mathbf{A(f_{x_i})}$ operators the Grover-like operator should be constructed using:\n", - "\n", - "$$\\mathbf{G}(\\mathbf{A(f_{x_i})}) = \\mathbf{A(f_{x_i})} \\left(\\hat{I} - 2|0\\rangle\\langle 0|\\right) \\mathbf{A(f_{x_i})}^{\\dagger}\\left(\\hat{I} - 2|\\Psi_0\\rangle\\langle \\Psi_0|\\right)$$" - ] - }, - { - "cell_type": "markdown", - "id": "94318e62", - "metadata": {}, - "source": [ - "In the case of our **QQuantLib** the Grover-like operator building can be done in an easy way using the **grover** function from *QQuantLib.AA.amplitude_amplification* module. The input of this function will be the following attribtues of the encoding class:\n", - "\n", - "* *oracle*\n", - "* *target*\n", - "* *index*" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "f42422c5", - "metadata": {}, - "outputs": [], - "source": [ - "from QQuantLib.AA.amplitude_amplification import grover" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b6d320c1", - "metadata": {}, - "outputs": [], - "source": [ - "operator_G = grover(\n", - " oracle = encoding_object.oracle,\n", - " target = encoding_object.target,\n", - " index = encoding_object.index\n", - ")" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "bac5d25a", - "metadata": {}, - "outputs": [], - "source": [ - "%qatdisplay operator_G --depth 2 --svg" - ] - }, - { - "cell_type": "markdown", - "id": "46e09369", - "metadata": {}, - "source": [ - "## 8. Amplitude Estimation Algorithm\n", - "\n", - "For each of the 3 $\\mathbf{A(f_{x_i})}$ operators and their correspondent Grover-like operators $\\mathbf{G}(\\mathbf{A(f_{x_i})})$ the desired *Amplitude Estimation* algorithm can be used. \n", - "\n", - "An amplitude estimation algorithm, usually, returns the probabiliy of getting the state $|\\Psi_0\\rangle$, so a post post-proccesing, for getting the estimator of the integral, should be done, as explained in section 6:\n", - "\n", - "$$\\tilde{S}_{[a,b]} = \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}\\right)$$\n", - "\n", - "\n", - "In case the algorithm returns some different value an additional post-processing should be done in order to get the $\\tilde{S}_{[a,b]}$ integral estimator." - ] - }, - { - "cell_type": "markdown", - "id": "794eee23", - "metadata": {}, - "source": [ - "In the case of **QQuantLib** steps 7 and 8 can be done straightforward using the different *Amplitude Estimation* classes of the modules from package *QQuantLib.AE*. \n", - "In fact using the **AE** class from *QQuantLib.AE.ae_class* we can acces to all the implemented algorithms in an easy way" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e84ed34a", - "metadata": {}, - "outputs": [], - "source": [ - "#This cell loads the QLM solver.\n", - "#QLMaaS == False -> uses PyLinalg\n", - "#QLMaaS == True -> try to use LinAlg (for using QPU as CESGA QLM one)\n", - "from QQuantLib.utils.qlm_solver import get_qpu\n", - "linalg_qpu = get_qpu('c')" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6db22fa8", - "metadata": {}, - "outputs": [], - "source": [ - "from QQuantLib.AE.ae_class import AE" - ] - }, - { - "cell_type": "markdown", - "id": "e16b209a", + "id": "ec6cc887-62d5-49db-95c7-95ef59220534", "metadata": {}, "source": [ - "First we create a base configuration dictionary for the **AE** object" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "57a5da80", - "metadata": {}, - "outputs": [], - "source": [ - "#AE base configuration dictionary\n", - "ae_dictionary = {\n", - " #QPU\n", - " 'qpu': linalg_qpu,\n", - " #Multi controlled decomposition\n", - " 'mcz_qlm': False, \n", - " \n", - " #shots\n", - " 'shots': None,\n", - " \n", - " #MLAE\n", - " 'schedule': [\n", - " [],\n", - " []\n", - " ],\n", - " 'delta' : None,\n", - " 'ns' : None,\n", - " \n", - " #CQPEAE\n", - " 'auxiliar_qbits_number': None,\n", - " #IQPEAE\n", - " 'cbits_number': None,\n", - " #IQAE & RQAQE\n", - " 'epsilon': None,\n", - " #IQAE\n", - " 'alpha': None,\n", - " #RQAE\n", - " 'gamma': None,\n", - " 'q': None\n", - "}" - ] - }, - { - "cell_type": "markdown", - "id": "349eda8e", - "metadata": {}, - "source": [ - "We are going to use the **IQAE** algorithm for solving the integral" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "73c20e30", - "metadata": {}, - "outputs": [], - "source": [ - "ae_dictionary.update({\n", - " 'ae_type': 'IQAE',\n", - " 'epsilon' : 0.001,\n", - " 'alpha': 0.05\n", - "})\n", - "ae_object = AE(\n", - " oracle=encoding_object.oracle,\n", - " target=encoding_object.target,\n", - " index=encoding_object.index,\n", - " **ae_dictionary\n", - ")\n", - "ae_object.run()\n", - "#We need to post-procces the return in order to get the correct estimator\n", - "ae_estimator_S = (b-a)*normalization*(2**n*np.sqrt(ae_object.ae_pdf))/2**n" - ] - }, - { - "cell_type": "markdown", - "id": "8e61b5e8", - "metadata": {}, - "source": [ - "In the case of the **IQAE** algorithm the lower and te upper limits for the amplitude is provided!!" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "6e1d389e", - "metadata": {}, - "outputs": [], - "source": [ - "ae_estimator_S" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "21698868", - "metadata": {}, - "outputs": [], - "source": [ - "print('Integral in first domain: {}'.format(sin_integral(a,b)))" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "a60e682a", - "metadata": {}, - "outputs": [], - "source": [ - "print(\"Is it the integral between the bound limits of the estimator: {}\".format(\n", - " (sin_integral(a,b) < ae_estimator_S['ae_u'])[0] & (sin_integral(a,b) > ae_estimator_S['ae_l'])[0]\n", - "))\n" - ] - }, - { - "cell_type": "markdown", - "id": "197ae2a6", - "metadata": {}, - "source": [ - "In the case of the **RQAE** algorithm the returned value is directly the **AMplitude** and not he probability (like in the ofhter methods implemented in the libary) so an additionall post-procces is mandatory. In this case:\n", - "\n", - "$$\\tilde{S}_{[a,b]} = \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\mathbf{P}[|\\Psi_0\\rangle]\\right)$$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "28107426", - "metadata": {}, - "outputs": [], - "source": [ - "ae_dictionary.update({\n", - " 'ae_type': 'RQAE',\n", - " 'epsilon' : 0.001,\n", - " 'alpha': None,\n", - " 'gamma': 0.05\n", - "})\n", - "ae_object = AE(\n", - " oracle=encoding_object.oracle,\n", - " target=encoding_object.target,\n", - " index=encoding_object.index,\n", - " **ae_dictionary\n", - ")\n", - "ae_object.run()\n", - "\n", - "#We need to post-procces the return in order to get the correct estimator.\n", - "#In the RQAE case the amplitude instead of the probability is returned!!!\n", - "ae_estimator_S = (b-a)*normalization*2**n*ae_object.ae_pdf/2**n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b29f2100", - "metadata": {}, - "outputs": [], - "source": [ - "print(\"Is it the integral between the bound limits of the estimator: {}\".format(\n", - " (sin_integral(a,b) < ae_estimator_S['ae_u'])[0] & (sin_integral(a,b) > ae_estimator_S['ae_l'])[0]\n", - "))\n" - ] - }, - { - "cell_type": "markdown", - "id": "f3891ada", - "metadata": {}, - "source": [ - "#### q_solve_integral\n", - "\n", - "The *q_solve_integral* function from *QQuantLib.finance.quantum_integration* module allows to solve integrals by *Amplitude Estimation* techniques of inputs arrays. \n", - "\n", - "The input of the function will be a complete dictionary with the configuration of the *amplitude estimation* algorithm and information about the arrays:\n", - "\n", - "* array_function : numpy array with the desired array function for Riemann sum\n", - "* array_probability : numpy array a probability distribtuion for computation of expected values. In our case this will be None.\n", - "* encoding : int for selecting the encoding. In the case of the benchmark will be 2.\n", - "\n", - "The outputs of the *q_solve_integral* function are:\n", - "* ae_estimation: pandas DataFrame with the desired integral computation and the upper and lower limits if applied.\n", - "* solver_ae: objet based on the AE class\n", - "\n", - "The *q_solve_integral* returns directly the integral of the input function directly. So the returned will be: $2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}$ (or the $2^{n} \\mathbf{P}[|\\Psi_0\\rangle]$ for the **RQAE** algorithm)\n", - "\n", - "\n", - "**BE AWARE!!**\n", - "\n", - "This is why we kept the $2^n$ in the integral computation, for taking advance of the implemented integral in the *q_solve_integral* function!!\n", - "\n", - "$$\n", - "\\tilde{S}_{[a,b]} = \\left\\{\n", - "\\begin{array}{ll}\n", - " \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}\\right) & For \\; probability \\; measurements \\\\\n", - " \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\mathbf{P}[|\\Psi_0\\rangle]\\right) & For \\; amplitude \\; measurements \\\\\n", - "\\end{array} \n", - "\\right.\n", - "$$" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "868865f1", - "metadata": {}, - "outputs": [], - "source": [ - "from QQuantLib.finance.quantum_integration import q_solve_integral" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ca9d8a30", - "metadata": {}, - "outputs": [], - "source": [ - "#AE base configuration dictionary\n", - "ae_dictionary = {\n", - " #QPU\n", - " 'qpu': linalg_qpu,\n", - " #Multi controlled decomposition\n", - " 'mcz_qlm': False, \n", - " \n", - " #shots\n", - " 'shots': None,\n", - " \n", - " #MLAE\n", - " 'schedule': [\n", - " [],\n", - " []\n", - " ],\n", - " 'delta' : None,\n", - " 'ns' : None,\n", - " \n", - " #CQPEAE\n", - " 'auxiliar_qbits_number': None,\n", - " #IQPEAE\n", - " 'cbits_number': None,\n", - " #IQAE & RQAQE\n", - " 'epsilon': None,\n", - " #IQAE\n", - " 'alpha': None,\n", - " #RQAE\n", - " 'gamma': None,\n", - " 'q': None\n", - "}\n", - "#IQAE configuration\n", - "ae_dictionary.update({\n", - " 'ae_type': 'IQAE',\n", - " 'epsilon' : 0.001,\n", - " 'alpha': 0.05\n", - "})\n", - "\n", - "encoding_dict = {\n", - " \"array_function\" : f_norm_x,\n", - " \"array_probability\" : None,\n", - " \"encoding\" : 2 \n", - "}\n", - "\n", - "ae_dictionary.update(encoding_dict)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "99792f3d", - "metadata": {}, - "outputs": [], - "source": [ - "solution, solver_object = q_solve_integral(**ae_dictionary)\n", - "#Integral of the input array is returned so only normalization has to be took into account\n", - "ae_estimator_S = normalization*solution*(b-a)/2**n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "e3f374ef", - "metadata": {}, - "outputs": [], - "source": [ - "ae_estimator_S" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "bebd04a9", - "metadata": {}, - "outputs": [], - "source": [ - "print(\"Is it the integral between the bound limits of the estimator: {}\".format(\n", - " (sin_integral(a,b) < ae_estimator_S['ae_u'])[0] & (sin_integral(a,b) > ae_estimator_S['ae_l'])[0]\n", - "))" - ] - }, - { - "cell_type": "markdown", - "id": "4ea9a2d5", - "metadata": {}, - "source": [ - "When using the **RQAE** algorithm *q_solve_integral* deals with the internal normalisations and the integral of the input array is returned in an transparent way." - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cad85df4", - "metadata": {}, - "outputs": [], - "source": [ - "#AE base configuration dictionary\n", - "ae_dictionary = {\n", - " #QPU\n", - " 'qpu': linalg_qpu,\n", - " #Multi controlled decomposition\n", - " 'mcz_qlm': False, \n", - " \n", - " #shots\n", - " 'shots': None,\n", - " \n", - " #MLAE\n", - " 'schedule': [\n", - " [],\n", - " []\n", - " ],\n", - " 'delta' : None,\n", - " 'ns' : None,\n", - " \n", - " #CQPEAE\n", - " 'auxiliar_qbits_number': None,\n", - " #IQPEAE\n", - " 'cbits_number': None,\n", - " #IQAE & RQAQE\n", - " 'epsilon': None,\n", - " #IQAE\n", - " 'alpha': None,\n", - " #RQAE\n", - " 'gamma': None,\n", - " 'q': None\n", - "}\n", - "#IQAE configuration\n", - "ae_dictionary.update({\n", - " 'ae_type': 'RQAE',\n", - " 'epsilon' : 0.001,\n", - " 'gamma': 0.05,\n", - " 'q' : 1.2\n", - "})\n", - "\n", - "encoding_dict = {\n", - " \"array_function\" : f_norm_x,\n", - " \"array_probability\" : None,\n", - " \"encoding\" : 2 \n", - "}\n", - "\n", - "ae_dictionary.update(encoding_dict)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "ac004a01", - "metadata": {}, - "outputs": [], - "source": [ - "solution, solver_object = q_solve_integral(**ae_dictionary)\n", - "#Integral of the input array is returned so only normalization has to be took into account\n", - "ae_estimator_S = normalization*solution*(b-a)/2**n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c778e480", - "metadata": {}, - "outputs": [], - "source": [ - "ae_estimator_S" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "c6be726d", - "metadata": {}, - "outputs": [], - "source": [ - "print(\"Is it the integral between the bound limits of the estimator: {}\".format(\n", - " (sin_integral(a,b) < ae_estimator_S['ae_u'])[0] & (sin_integral(a,b) > ae_estimator_S['ae_l'])[0]\n", - "))" - ] - }, - { - "cell_type": "markdown", - "id": "9c1a0dd0", - "metadata": {}, - "source": [ - "## 9. Getting the metrics\n", - "\n", - "Once the amplitude estimator of the integral, $\\tilde{S}$, is obtained following metrics should be computed:\n", - "\n", - "* *Absolute error* between the ae estimator and the exact integral to compute: $\\epsilon = |\\tilde{S} - (\\cos(a)-\\cos(b))|$\n", - "\n", - "* *Oracle calls*: total number of calls of the operator $\\mathbf{A}$ (the number of shots should be taking into account in this calculation).\n", - "\n", - "* *Elapsed time*:The complete time for getting $\\tilde{S}$. This time will include all the complete steps this is (it is not necesary have times of each individual step):\n", - " * Discretization time\n", - " * Creation of operator $\\mathbf{A}$\n", - " * Creation of the Grover-like operator: $\\mathbf{G}$\n", - " * Execution of the *amplitude Estimation* algorithm\n", - " * Post-processing execution\n", - " * Metrics calculation." - ] - }, - { - "cell_type": "markdown", - "id": "9e7a56cb", - "metadata": {}, - "source": [ - "\n", - "The *sine_integral* function from *AmplitudeEstimation/ae_sine_integral* allows the user perform one complete computation of the sine integral by **AE** techniques. The inputs are:\n", - "\n", - "* *n_qbits* : number of qubits for interval discretization.\n", - "* *interval* : int for selecting the integration interval.\n", - "* *ae_dictionary* : python dictionary with the complete configuration of the *Amplitude Estimation* algorithm.\n", - "\n", - "The return is a panda DataFrames with the complete information of the benchamrk.\n" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3a802720", - "metadata": {}, - "outputs": [], - "source": [ - "#AE base configuration dictionary\n", - "ae_dictionary = {\n", - " #QPU\n", - " 'qpu': 'c',\n", - " #Multi controlled decomposition\n", - " 'mcz_qlm': False, \n", - " \n", - " #shots\n", - " 'shots': None,\n", - " \n", - " #MLAE\n", - " 'schedule': [\n", - " [],\n", - " []\n", - " ],\n", - " 'delta' : None,\n", - " 'ns' : None,\n", - " \n", - " #CQPEAE\n", - " 'auxiliar_qbits_number': None,\n", - " #IQPEAE\n", - " 'cbits_number': None,\n", - " #IQAE & RQAQE\n", - " 'epsilon': None,\n", - " #IQAE\n", - " 'alpha': None,\n", - " #RQAE\n", - " 'gamma': None,\n", - " 'q': None\n", - "}\n", - "#IQAE configuration\n", - "ae_dictionary.update({\n", - " 'ae_type': 'IQAE',\n", - " 'epsilon' : 0.001,\n", - " 'alpha': 0.05\n", - "})" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "2b30ff7e", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "from ae_sine_integral import sine_integral" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "170d2a4b", - "metadata": {}, - "outputs": [], - "source": [ - "#first interval integral benchmarking\n", - "pdf = sine_integral(4, 0, ae_dictionary)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cb94f6da", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "pdf[['integral_ae', 'integral_ae_l', 'integral_ae_u', 'exact_integral', 'absolute_error_sum']]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cf95c732", - "metadata": {}, - "outputs": [], - "source": [ - "#second integral is negative\n", - "pdf = sine_integral(4, 1, ae_dictionary)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "cd7df328", - "metadata": { - "scrolled": true - }, - "outputs": [], - "source": [ - "pdf[['integral_ae', 'integral_ae_l', 'integral_ae_u', 'exact_integral', 'absolute_error_sum']]" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "3afd06ee", - "metadata": {}, - "outputs": [], - "source": [ - "# Only RQAE can deal with negative integrals\n", - "#RQAE configuration\n", - "ae_dictionary.update({\n", - " 'ae_type': 'RQAE',\n", - " 'epsilon' : 0.001,\n", - " 'gamma': 0.05,\n", - " 'q' : 2.0,\n", - " 'alpha': None\n", - "})\n", - "#second integral is negative\n", - "pdf = sine_integral(4, 1, ae_dictionary)" - ] - }, - { - "cell_type": "code", - "execution_count": null, - "id": "b9be8c50", - "metadata": {}, - "outputs": [], - "source": [ - "pdf[['integral_ae', 'integral_ae_l', 'integral_ae_u', 'exact_integral', 'absolute_error_sum']]" - ] - }, - { - "cell_type": "markdown", - "id": "d85f36ed", - "metadata": {}, - "source": [ - "## 10. Benchmark Summary\n", - "\n", - "We have all the pieces needed for do a complete benchmarking. Now we summarize the procedure:\n", - "\n", - "1. Execute following steps for each one of the 3 integration described intervals in Sub-Section 6.1\n", - " 1. Execute the following steps from n = 4 to the maximum nubmer of qubits. For each n execute following steps 10 times.\n", - " 1. Create the domain discretization (sub-section 6.2)\n", - " 2. Create the array with the correspondient sine function discretization (sub-section 6.3)\n", - " 3. Compute normalization of the array (sub-section 6.4)\n", - " 4. Create the $\\mathbf{A}$ oracle operator for encoding the array (see sub-section 6.5)\n", - " 5. Create the correspondient Grover-like operator, $\\mathbf{G}(\\mathbf{A(f_{x_i})})$, from $\\mathbf{A}$ (see section 7)\n", - " 6. Execute the tested *amplitude estimation* algorithm properly confured using the operators $\\mathbf{G}(\\mathbf{A(f_{x_i})})$, and $\\mathbf{A}$ (see section 8)\n", - " 7. With the *amplitude estimation* algorithm result compute the integral (see section 8)\n", - " 8. Compute the desired metrics (see section 9).\n", - " 2. For each of the n qbits tested computed the mean, standard deviation, minimum and maximum of the 10 values of each metric.\n", - "2. For each one of the interval integration return the computed metric statistics for each number of qbits. \n", - "\n", + "## 6. About the AE algorithms.\n", "\n", - "If a configuration parameter of the *amplitude estimation* algorithm want to be benchmarked for each posible values the complete above procedure should be executed.\n", + "Different **AE** algorithms can be used for solving the **AE** kernel. In the present library, the following ones were implemented:\n", "\n", - "Generating a report using the different results of the benchmarking." + "* Classical Quantum Phase Estimation (**CQPEAE**): based on *Brassard et al (2022). Quantum Amplitude Amplification and Estimation*. This is the presented method in section 4. **Only estimates positives $a$**.\n", + "* Iterative Quantum Phase Estimation (**IQPEAE**): based on *Kitaev (1995) Quantum measurements and the Abelian Stabilizer Problem*. The method is similar to **CQPEAE** but the **QFT** is performed using a single qubit. This method generates a circuit with a lower number of qubits but with deeper ones. **Only estimates positives $a$**\n", + "* Maximum Likelihood Amplitude Estimation (**MLAE**): based on *Suzuki et al (2020) Amplitude estimation without phase estimation*. In this case direct quantum circuits with different powers of the corresponding Grover operator $\\mathbf{G}(\\mathbf{A})$ (see section 3) are executed and measured. The different measurements are post-processed using maximum likelihood techniques for getting an estimation of $a$. **Only estimates positives $a$**\n", + "* Iterative Quantum Amplitude Estimation (**IQAE**): based on *Grinko et al (2021) Iterative quantum amplitude estimation*. This algorithm uses the Grover operator $\\mathbf{G}(\\mathbf{A})$ in an iterative way. In each step a differing number of applications of $\\mathbf{G}(\\mathbf{A})$ is computed based on the results obtained in the step before. This algorithm returns a $\\alpha$ level confidence interval for the estimation of $a$. So algorithm provides a minimum and maximum of $a$. The length of the interval should be fixed a priori. For lower lengths, the algorithm needs more steps and deeper circuits. **Only estimates positives $a$**\n", + "* Real Quantum Amplitude Estimation (**RQAE**): based on *Manzano et al (2023) Real Quantum Amplitude Estimation*. An iterative algorithm where in each step the number of applications of the Grover operator $\\mathbf{G}(\\mathbf{A})$ is computed using the results from before steps. This algorithm returns a $\\alpha$ level confidence interval for the estimation of $a$. So algorithm provides a minimum and maximum of $a$. The length of the interval should be fixed a priori. For lower lengths, the algorithm needs more steps and deeper circuits. **This algorithm allows us to estimate non-positive $a$**." ] } ], @@ -1369,7 +171,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.9.9" + "version": "3.11.7" } }, "nbformat": 4, diff --git a/tnbs/BTC_02_AE/QQuantLib/notebooks/02_AmplitudeEstimationBTC.ipynb b/tnbs/BTC_02_AE/QQuantLib/notebooks/02_AmplitudeEstimationBTC.ipynb new file mode 100644 index 0000000..8c186b5 --- /dev/null +++ b/tnbs/BTC_02_AE/QQuantLib/notebooks/02_AmplitudeEstimationBTC.ipynb @@ -0,0 +1,1317 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "0618c8ae", + "metadata": {}, + "source": [ + "# The Amplitude Estimation Benchmark Test Case\n", + "\n", + "In the **01_AmplitudeEstimationKernel.ipynb** notebook the **AE kernel** was explained. In this notebook the associated benchmark test case (**BTC**) is presented. " + ] + }, + { + "cell_type": "code", + "execution_count": 1, + "id": "7071be0e", + "metadata": {}, + "outputs": [], + "source": [ + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "%matplotlib inline" + ] + }, + { + "cell_type": "markdown", + "id": "011635f6", + "metadata": {}, + "source": [ + "## 1. Description of the problem\n", + "\n", + "The benchmark problem is the computation of the integral of a function $f(x)$ in a closed interval $[a,b] \\subset \\mathbf{R}$:\n", + "\n", + "$$I = \\int_a^bf(x)dx$$\n", + "\n", + "In particular we propose to use $f(x) = \\sin x$, whose integral is very well known:\n", + "\n", + "$$I = \\int_a^{b}\\sin(x)dx = -\\cos x |_a^b = \\cos(a)-\\cos(b)$$\n", + "\n", + "Additionally, 2 different integration intervals will be used:\n", + "\n", + "1. $[0, \\frac{3\\pi}{8}]$ : $I = \\int_0^{\\frac{3\\pi}{8}}\\sin(x)dx = 0.6173165676349102$\n", + "2. $[\\pi, \\frac{5\\pi}{4}]$ : $I = \\int_{\\pi}^{\\frac{5\\pi}{4}}\\sin(x)dx = -0.2928932188134523$\n" + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "871a0510", + "metadata": {}, + "outputs": [], + "source": [ + "def sin_integral(a,b):\n", + " return np.cos(a)-np.cos(b)" + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "c7aec8ad", + "metadata": {}, + "outputs": [], + "source": [ + "start = [0.0, 3.0*np.pi/4.0, np.pi]\n", + "end = [3.0*np.pi/8.0, 9.0*np.pi/8.0, 5.0*np.pi/4.0]" + ] + }, + { + "cell_type": "markdown", + "id": "148bbae2", + "metadata": {}, + "source": [ + "#### $1^{st}$ Domain: $[0, \\frac{3\\pi}{8}]$ \n", + "\n", + "In this domain the function is strictly positive defined." + ] + }, + { + "cell_type": "code", + "execution_count": 4, + "id": "29114c8f", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 4, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAAA97klEQVR4nO3deVxVdeLG8eeyCwqKC6Ig4r4vgJA6tkdZo9mmpbmlFdNqtunYZDnNWDY106YtbmlmpFnTlFnMVOaWC+KK+wayiGCyiJfl3vP7g4nfmFpcBM7l3s/79bp/dDyH+9xveO7jWb7HYhiGIQAAAJN4mB0AAAC4N8oIAAAwFWUEAACYijICAABMRRkBAACmoowAAABTUUYAAICpKCMAAMBUXmYHqAq73a7MzEw1atRIFovF7DgAAKAKDMNQYWGhWrVqJQ+Pix//qBdlJDMzU+Hh4WbHAAAA1ZCenq6wsLCL/nm9KCONGjWSVPFhAgMDTU4DAACqoqCgQOHh4ZXf4xdTL8rIz6dmAgMDKSMAANQzv3WJBRewAgAAU1FGAACAqSgjAADAVJQRAABgKsoIAAAwFWUEAACYijICAABMRRkBAACmoowAAABTUUYAAICpKCMAAMBUlBEAAGAqyggAAG7KWmbTv7ZnasLCzTp1ptS0HPXiqb0AAKBmGIahLcd+0oqtx/XFjiwVWsslSf/anqmxA9qakokyAgCAG8g4fVYrko9r+dbjOpZXXLm8deMGujWqta7s3Ny0bJQRAABclLXMpq93Z2t58nGtPZgrw6hYHuDjqcE9Q3VbVJjiIoPl4WExNSdlBAAAF7M7M18fb07XpykZKvjvaRhJ6t+uqe6ICdMNPVrK38d5KoDzJAEAANVWYC3T59sylbg5XTsz8iuXt27cQLdHh+n26DCFB/ubmPDiKCMAANRThmFo+/F8Ld2Yps+3Z+psmU2S5OPpoeu6h+jOfuEa2L6Z6adhfgtlBACAeqaopFyfpmRo6cY0pWYVVC7v0KKh7uwXrlujwhQc4GNiQsdQRgAAqCf2ZBXogx+P6bOUDJ0p/e9REC8P/b5nqO6Ka6OYiCayWJz7KMiFUEYAAHBipeV2fbUrS4s2HFPysZ8ql7dvHqCRcRG6Laq1GvvXn6MgF0IZAQDACWXln9WHG9O0dFOacosqZkf18rDo+u4tNeqyNurfrmm9PApyIZQRAACchGEY2nz0Jy1cf0Rf7z4hm71iYpCQQF+NjI3QXbHhahHoZ3LKmkcZAQDAZNYymz7fnqmF646ec0FqXGSwxg5oq+u6hcjb03UfJ0cZAQDAJDkFVi3+8ZiWbEyrfFCdn7eHbunbWmMHtFWXloEmJ6wblBEAAOrYrox8zV97RP/akakyW8WpmNaNG2h0/wiNiAlXk3p0W25NoIwAAFAH7HZD3+7N0XtrDmvjkVOVy2Mimuie30UqvluIvFz4VMyvoYwAAFCLrGU2fbL1uOatPaLDJ89Iqrgr5qZeobpnYKR6hzc2N6AToIwAAFALTp0p1aINR7Vow7HK60Ea+XlpZFwbjRvQVqFBDUxO6DwoIwAA1KD0U8Wau+awEreky1pmlySFNWmgewZGani/cDX05av3lxgRAABqwO7MfL29+rC+3JGp/04Pop6tg3Tf5e00uEdLt70epCooIwAAVJNhGNp45JTmfH9Iq/efrFx+eafmSri8nfq3d51ZUmsTZQQAAAfZ7Yb+szdHs78/qJS005IkD4s0pHcr3X95e3Vr5R7zg9QUyggAAFVksxv6Ykem5nx/SHuzCyVVPDV3eEyY7hvUXm2a+pucsH6ijAAA8BtKy+36NOW45nx/SEfziiVJDX29NOqyNprwu0i1aOR6z4upS5QRAAAuoqTcpo+3HNfb3x9SxumzkqTG/t66Z2CkxvZvqyB/b5MTugbKCAAAv2Ats2nppjS9vfqQThSUSJKaN/LVfYPaaWRcGwVwe26NYjQBAPgva5lNSzZWlJCThRUlpGWgnxKuaKc7Y9vIz9vT5ISuiTICAHB7Z0ttWrLxmN5efVi5RRUlpHXjBnrgqva6PTpMvl6UkNpEGQEAuK2fj4TM+f7QOSXk4as76NaoMPl4MVFZXaCMAADcTkm5TYmb0/XWdwcrrwmhhJiHMgIAcBtlNruWJx/XG/85oMx8qySpVZCfHr6mo26jhJiGMgIAcHk2u6F/bsvQP/59QGmnKuYJCQn01UNXddDwfuFcE2IyyggAwGXZ7YZW7c7Wq0n7dTCnSJLUrKGPHriyg0bGcXeMs6CMAABcjmEY+uFArl7+eq92ZRRIkoIaeOv+K9pp3IC28vfh68+Z8H8DAOBStqb9pFmr9urHw6ckSQE+npowqJ0mDopUoB8zpjojyggAwCXsP1Gol7/ep6TUE5IqHmA3+rIIPXBlezVt6GtyOvwayggAoF7LPH1Wf0/ar0+2HpfdkDws0h3R4Xr02o5q1biB2fFQBZQRAEC9dLq4VLO/P6SF64+qtNwuSRrco6Uej++sDi0ampwOjqCMAADqFWuZTe+vP6q3vjuoAmu5JCkuMlhPD+6iqDZNTE6H6qCMAADqBbvd0GfbMvTKN/uVcfqsJKlLy0Z6enAXXdmpuSwWi8kJUV2UEQCA01t3MFd/+XKPUrMqbtMNDfLT4/GddUvf1vL0oITUd5QRAIDT2n+iUDNX7tF3+05Kkhr5eemBKzto/MC2TFjmQqo1Cf/s2bMVGRkpPz8/RUdHa82aNb+6/pIlS9S7d2/5+/srNDRU48ePV15eXrUCAwBc38nCEk1dsVM3/OMHfbfvpLw8LBo3oK1WP3mV/nBle4qIi3G4jCQmJmrSpEmaNm2aUlJSNGjQIA0ePFhpaWkXXH/t2rUaM2aMJkyYoN27d2vZsmXavHmzJk6ceMnhAQCuxVpm01vfHdSVL3+npZvSZDekG7q3VNLkK/Tc0O4KDvAxOyJqgcUwDMORDeLi4hQVFaU5c+ZULuvatauGDRummTNnnrf+3/72N82ZM0eHDh2qXPbGG29o1qxZSk9Pr9J7FhQUKCgoSPn5+QoMDHQkLgCgHjAMQ//akaWXvtpbeXFq77AgPfP7burXNtjkdKiuqn5/O3RkpLS0VMnJyYqPjz9neXx8vNavX3/BbQYMGKDjx49r5cqVMgxDJ06c0PLly3XTTTc58tYAABe1Lf20bpuzXo8sTVHG6bMKDfLTP0b00acPDKSIuAmHLmDNzc2VzWZTSEjIOctDQkKUnZ19wW0GDBigJUuWaMSIEbJarSovL9fQoUP1xhtvXPR9SkpKVFJSUvnfBQUFjsQEANQD2flWzVq1VytSMiRJ/j6e+sMV7TVxUDs18OGaEHdSrQtYf3kvt2EYF72/OzU1VY888oieffZZJScna9WqVTpy5IgSEhIu+vNnzpypoKCgyld4eHh1YgIAnJC1zKbX/3NAV/3t+8oicnt0mL574ko9fE1HiogbcuiakdLSUvn7+2vZsmW65ZZbKpc/+uij2rZtm1avXn3eNqNHj5bVatWyZcsql61du1aDBg1SZmamQkNDz9vmQkdGwsPDuWYEAOoxwzC0cme2/rpyT+V1IdERTTR9SDf1CmtsbjjUiqpeM+LQaRofHx9FR0crKSnpnDKSlJSkm2+++YLbFBcXy8vr3Lfx9KxovRfrQb6+vvL15QmLAOAq9mQV6Pl/7daPh09Jqpi0bOqNXTWkVygzp8LxSc8mT56s0aNHKyYmRv3799e7776rtLS0ytMuU6dOVUZGhhYtWiRJGjJkiO69917NmTNH119/vbKysjRp0iTFxsaqVatWNftpAABO5aczpXo1ab+WbDwmuyH5enno/iva6w9XtOd0DCo5XEZGjBihvLw8zZgxQ1lZWerRo4dWrlypiIgISVJWVtY5c46MGzdOhYWFevPNN/X444+rcePGuvrqq/XSSy/V3KcAADgVm93QR5vT9PLX+3S6uEySdGPPlpo6uKvCg/1NTgdn4/A8I2ZgnhEAqD+Sj53S9M93a1dGxZ2QnUMaafrQbhrQvpnJyVDXauWaEQAALuZkYYle/GqvPtl6XFLFc2QmX9dJoy+LkJdntW7ehJugjAAALkm5za4lG9P0t2/2qdBaLkkaEROuJ2/orGYNuRkBv40yAgCotuRjp/Snz3YrNavilEzP1kGacXN39W3TxORkqE8oIwAAh+UVVZySWZZccUomqIG3nry+s+6KbSNPD27VhWMoIwCAKrPbDX20OV0vrdqr/LMVd8mMiAnXUzd0VlNOyaCaKCMAgCrZlZGvaZ/t0vb005KkrqGBemFYD0VHcEoGl4YyAgD4VYXWMr3yzX4t2nBUdkNq6Oulx+O5SwY1hzICALggwzD01a5sPf+v3TpRUPG8sCG9W+lPN3VVi0A/k9PBlVBGAADnST9VrD/9c5e+33dSktS2qb/+PKyHBnVsbnIyuCLKCACgUpnNrrlrjui1/+yXtcwuH08P/eHK9vrDle3l582zZFA7KCMAAEnS1rSf9McVO7U3u1CS1L9dU71wSw+1b97Q5GRwdZQRAHBzBdYyzVq1V0s2pskwpCb+3nrmpm66Naq1LBbmDEHto4wAgBtbtStb0z/fVXmB6u3RYfrjjV0VHOBjcjK4E8oIALih7Hyrnv3nLn2TekJSxQWqf72lpwZ04Mm6qHuUEQBwI3a7oSWb0jTrq70qLCmXl4dFCVe010NXd+ACVZiGMgIAbuJgTpGmrtihzUd/kiT1CW+sF2/rqS4tA01OBndHGQEAF1dms+ud1Yf0+n8OqtRml7+Pp566vrNG92/LQ+3gFCgjAODCdhw/raeW76i8XffKzs31l1t6qnXjBiYnA/4fZQQAXJC1zKa/J+3Xe2sOy/7f23WnD+mum/u04nZdOB3KCAC4mI2H8/T0Jzt0NK9YkjS0dytNH9JNTRv6mpwMuDDKCAC4iKKScr341R598GOaJCkk0Fd/GdZT13YLMTkZ8OsoIwDgAtYcOKkpn+xUxumzkqS7YsM19cauCvTzNjkZ8NsoIwBQjxVYy/TXL/foo83pkqSwJg300m29NJDJy1CPUEYAoJ76fl+Opq7Yqax8qyRpbP8IPXVDFwX4smtH/cJvLADUMwXWMr3wRao+3nJckhTR1F+zbuuluHZNTU4GVA9lBADqkdX7T2rKJzuUlW+VxSKNHxCpJ6/vrAY+TOWO+osyAgD1QKG1TC98sUeJWyquDWnb1F+zbu+t2Mhgk5MBl44yAgBObt3BXD21fIcyTp+VxSKNG9BWT13fhaMhcBmUEQBwUmdKyvXiV3u1+MdjkqQ2wf56+XauDYHroYwAgBPafPSUHv94u9JOVcyiOvqyCE0ZzJ0ycE38VgOAE7GW2fTqf58pYxhSqyA/zbq9t37XkXlD4LooIwDgJHZl5Gvyx9u0/0SRJOmO6DD9aUg3ZlGFy6OMAIDJymx2zf7ukN749oDK7YaaNfTVi7fyTBm4D8oIAJjo0MkiTU7cpu3H8yVJN/ZsqReG9VRwgI/JyYC6QxkBABPY7YYW/3hMM7/aI2uZXYF+Xppxcw/d3KeVLBaL2fGAOkUZAYA6lp1v1ZPLt2vNgVxJ0u86NNPLd/RSaFADk5MB5qCMAEAd+tf2TD3z2S7lny2Tr5eHpg7uojH928rDg6MhcF+UEQCoA/lnyzT9n7v02bZMSVLP1kH6+4g+6tCiocnJAPNRRgCglm04lKfHP96mzHyrPCzSg1d10CPXdJS3p4fZ0QCnQBkBgFpSWm7XK0n79O4PFROYtQn2199H9FZ0BA+3A/4XZQQAasHBnEI9+tE27c4skCQNjwnTs0O6qyHTuQPn4W8FANQgwzD0wY/H9MKXe1RSbldjf2+9eGsv3dCjpdnRAKdFGQGAGpJbVKKnlu/Qt3tzJEmDOjbT3+7orZBAP5OTAc6NMgIANeD7fTl6YtkO5RaVyMfLQ1Nu6KJxA7hlF6gKyggAXAJrmU0vrdqrBeuOSpI6hTTU63f1VZeWgeYGA+oRyggAVNP+E4V6ZGmK9mYXSpLGDWirKYO7yM/b0+RkQP1CGQEABxmGoSUb0/TnL1JVUm5Xs4Y+evn23rqqSwuzowH1EmUEABzw05lSPf3JDn2TekKSdHmn5vrbHb3UohEXqQLVRRkBgCracChPjyVuU3aBVd6eFj19QxfdMzCSi1SBS0QZAYDfUG6z67X/HNCb3x2UYUjtmgfo9Tv7qkfrILOjAS6BMgIAv+L4T8V69KNtSj72k6SKmVSfG9pd/j7sPoGawt8mALiIlTuzNOWTHSqwlquRr5f+emtPDendyuxYgMuhjADAL1jLbJrxRao+3JgmSeoT3lhv3NVX4cH+JicDXBNlBAD+x4EThXrowxTtO1Eoi0VKuKK9Jl/XSd6eHmZHA1wWZQQAVDF3yLItx/Xs57tkLauYO+TvI/poUMfmZkcDXB5lBIDbKyop17RPd+qf2zIlVTzg7pXhvZk7BKgjlBEAbm1XRr4e+nCrjuYVy9PDosfjOynh8vbMHQLUIcoIALdkGIYW/3hML3yxR6U2u1o3bqDX7+qj6Ihgs6MBbocyAsDt5J8t05RPduirXdmSpGu7huhvd/RSY38fk5MB7okyAsCtbE8/rYeWblX6qbPy9rRo6uCuGj+wrSwWTssAZqGMAHALhmHo/fVH9ZeVe1RmMxQe3EBv3hWl3uGNzY4GuD3KCACXl3+2TE8v36FVuytOy9zQvaVeur2Xghp4m5wMgCRVaxaf2bNnKzIyUn5+foqOjtaaNWt+df2SkhJNmzZNERER8vX1Vfv27TV//vxqBQYAR+w8nq/fv7FGq3Zny9vToueGdNOcu6MoIoATcfjISGJioiZNmqTZs2dr4MCBeueddzR48GClpqaqTZs2F9xm+PDhOnHihObNm6cOHTooJydH5eXllxweAC7ml3fLhDVpoLdGcloGcEYWwzAMRzaIi4tTVFSU5syZU7msa9euGjZsmGbOnHne+qtWrdKdd96pw4cPKzi4erfMFRQUKCgoSPn5+QoMDKzWzwDgPgqtZZqyYqe+3JElSYrvFqKX7+jN0RCgjlX1+9uh0zSlpaVKTk5WfHz8Ocvj4+O1fv36C27z+eefKyYmRrNmzVLr1q3VqVMnPfHEEzp79uxF36ekpEQFBQXnvACgKvZkFWjom+v05Y4seXlY9MxNXfXO6GiKCODEHDpNk5ubK5vNppCQkHOWh4SEKDs7+4LbHD58WGvXrpWfn58+/fRT5ebm6oEHHtCpU6cuet3IzJkz9fzzzzsSDQD08ZZ0/emzXSoptys0yE9vjoxSdEQTs2MB+A3VuoD1l/fjG4Zx0Xv07Xa7LBaLlixZotjYWN1444169dVXtXDhwoseHZk6dary8/MrX+np6dWJCcBNWMtsemr5dj21fIdKyu26olNzffnIIIoIUE84dGSkWbNm8vT0PO8oSE5OznlHS34WGhqq1q1bKygoqHJZ165dZRiGjh8/ro4dO563ja+vr3x9fR2JBsBNHck9oweWbNWerAJ5WKTJ13XSA1d24NkyQD3i0JERHx8fRUdHKykp6ZzlSUlJGjBgwAW3GThwoDIzM1VUVFS5bP/+/fLw8FBYWFg1IgNAhVW7sjT0jbXak1WgZg19tHhCnB66uiNFBKhnHD5NM3nyZM2dO1fz58/Xnj179NhjjyktLU0JCQmSKk6xjBkzpnL9kSNHqmnTpho/frxSU1P1ww8/6Mknn9Q999yjBg0a1NwnAeA2ymx2/eXLVCV8sFWFJeXq17aJvnxkkAZ2aGZ2NADV4PA8IyNGjFBeXp5mzJihrKws9ejRQytXrlRERIQkKSsrS2lpaZXrN2zYUElJSXr44YcVExOjpk2bavjw4XrhhRdq7lMAcBs5BVY99GGKNh09JUm67/J2evL6zvL2rNYlcACcgMPzjJiBeUYASNKPh/P00Icpyi0qUUNfL/3tjl66oUeo2bEAXERVv795Ng0Ap2cYht794bBmfb1PNruhLi0bafaoKLVr3tDsaABqAGUEgFMrtJbpyWX//5C7W/u21l9u6akGPp4mJwNQUygjAJzW/hOFSlicrMO5Z+TtadH0Id01Kq7NRec1AlA/UUYAOKXPt2fq6eU7dLbMptAgP80eFaW+bZjEDHBFlBEATqXMZtfMlXs1f90RSdLADk31+p191bQhEyECrooyAsBp5BRa9dCS/79t94Er2+vx+M7yZBIzwKVRRgA4heRjp/SHD7Yqp/Dn23Z764YeLc2OBaAOUEYAmMowDC3+8Zhm/CtV5XZDHVs01Dujo7ltF3AjlBEAprGW2TTt0136ZOtxSdJNvUI167ZeCvBl1wS4E/7GAzBF+qliJXyQrN2ZFU/bnTq4qyYOiuS2XcANUUYA1Lm1B3L18NKt+qm4TMEBPnrzrr4awEPuALdFGQFQZwzD0Ds/HNasVXtlN6TeYUGafXe0WjfmCd6AO6OMAKgTxaXlenL5Dn25I0uSdEd0mP48rIf8vJnWHXB3lBEAte5o7hndvzhZ+04UytvTomeHdNfdTOsO4L8oIwBq1ff7cvTI0hQVWMvVvJGv5oyKUkzbYLNjAXAilBEAtcIwDM1ZfUgvf71PhiH1bdNYb98drZBAP7OjAXAylBEANe6X14fcFRuu54Z2l68X14cAOB9lBECNSssr1n2Lt2hvdsX1Ic8N7a5RcRFmxwLgxCgjAGrM2gO5emjpVp0uLlOzhr56+26uDwHw2ygjAC6ZYRiat/aI/rpyT8X8IeGN9c7d0WoZxPUhAH4bZQTAJbGW2fTHT3dqxdYMSdJtUWH6yy3MHwKg6igjAKotK/+sEhYna/vxfHl6WDTtxq4aP7At84cAcAhlBEC1JB87pfsXb1VuUYka+3tr9sgoni8DoFooIwAc9vHmdD3z2S6V2uzq0rKR3hsTo/Bgf7NjAainKCMAqqzMZtdfvtyjheuPSpJu6N5SrwzvrQBfdiUAqo89CIAq+elMqR5aulXrDuZJkh67tpMevrqDPDy4PgTApaGMAPhN+08U6t5FW3Qsr1j+Pp56dXgf3dCjpdmxALgIygiAX/Xv1BOalLhNRSXlCmvSQHPHxqhLy0CzYwFwIZQRABf0ywfdxUUGa87d0QoO8DE7GgAXQxkBcB5rmU1Pf7JD/9yWKUkafVmEnh3STd6eHiYnA+CKKCMAznGiwKr7Fm3R9uP58vKoeNDd3ZfxoDsAtYcyAqDS9vTTum/xFp0oqJjIbM6oaPVv39TsWABcHGUEgCTpn9sy9NTyHSopt6tTSEPNHdNPbZoykRmA2kcZAdyc3W7o1aT9evO7g5Kka7q00D/u7KNGft4mJwPgLigjgBs7U1KuyR9v09e7T0iS7r+8nZ66oYs8mcgMQB2ijABuKuP0WU18f4v2ZBXIx9NDf721p26PDjM7FgA3RBkB3NDWtJ9036Jk5RaVqFlDH70zOlrREcFmxwLgpigjgJv557YMPbl8h0rLK564O3dsjMKacKEqAPNQRgA3Ybcb+vu/9+uNbysuVL22a4heu7MPT9wFYDr2QoAbKC4t1+Mfb9dXu7IlSfdf0U5PXc+FqgCcA2UEcHHZ+Vbdu2iLdmbky9vTopm39uJCVQBOhTICuLCdx/M1cdFmnSgoUXBAxYWq/dpyoSoA50IZAVzUql1ZmpS4TdYyuzq2aKh5Y5lRFYBzoowALsYwDM3+/pBe/nqfJOnyTs315si+CmRGVQBOijICuJCScpv+uGKXPtl6XJI0bkBbPXNTV3l5epicDAAujjICuIhTZ0qVsDhZm46ekqeHRc8N6abR/duaHQsAfhNlBHABh04W6Z6Fm3Usr1iNfL305qgoXdGpudmxAKBKKCNAPbf+YK4SPkhWgbVcYU0aaP64fuoU0sjsWABQZZQRoB5L3JymaZ/uUrndUFSbxnp3TIyaNfQ1OxYAOIQyAtRDdruhl77eq3dWH5YkDe3dSrNu7yU/b0+TkwGA4ygjQD1zttSmxxK3adXuiqndH72moyZd21EWC1O7A6ifKCNAPZJTYNXERVu043i+fDw9NOv2XhrWt7XZsQDgklBGgHpib3aB7lmwWZn5VjXx99a7Y2KY2h2AS6CMAPXA9/ty9NCHKSoqKVe75gFaMK6fIpoGmB0LAGoEZQRwch/8eEzTP98tm93QZe2C9c7dMQryZ2p3AK6DMgI4Kbvd0Myv9ui9NUckSbdFhWnmrT3l48XU7gBcC2UEcEJnS22alJiir3efkCQ9Ed9JD17VgTtmALgkygjgZE4Wlmjioi3ann5aPp4eevmOXrq5D3fMAHBdlBHAiRw4UahxCzYr4/RZ7pgB4DYoI4CTWPffZ8wUWssV2SxA88f1U2Qz7pgB4PooI4ATWLYlXVNX7FS53VC/tk307ugYNQnwMTsWANQJyghgIsMw9Pek/Xr924OSpCG9W+llnjEDwM1U6x7B2bNnKzIyUn5+foqOjtaaNWuqtN26devk5eWlPn36VOdtAZdSUm7T5I+3VxaRB69qr9dG9KGIAHA7DpeRxMRETZo0SdOmTVNKSooGDRqkwYMHKy0t7Ve3y8/P15gxY3TNNddUOyzgKvKLyzRm3iZ9mpIhTw+LXry1p568vos8PLh1F4D7sRiGYTiyQVxcnKKiojRnzpzKZV27dtWwYcM0c+bMi2535513qmPHjvL09NRnn32mbdu2Vfk9CwoKFBQUpPz8fAUGBjoSF3A66aeKNW7BJh06eUYNfb00e1SULu/U3OxYAFDjqvr97dCRkdLSUiUnJys+Pv6c5fHx8Vq/fv1Ft1uwYIEOHTqk6dOnV+l9SkpKVFBQcM4LcAXb00/rltnrdOjkGbUM9NOyhP4UEQBuz6EykpubK5vNppCQkHOWh4SEKDs7+4LbHDhwQFOmTNGSJUvk5VW162VnzpypoKCgyld4eLgjMQGnlJR6Qne++6Nyi0rVNTRQnz04UF1DOdIHANW6gPWXU1IbhnHBaaptNptGjhyp559/Xp06daryz586dary8/MrX+np6dWJCTiNRRuO6v7FW3S2zKbLOzXXx/dfppZBfmbHAgCn4NCtvc2aNZOnp+d5R0FycnLOO1oiSYWFhdqyZYtSUlL00EMPSZLsdrsMw5CXl5e++eYbXX311edt5+vrK19fX0eiAU7Jbjf04qq9eveHw5KkO/uF68/Desjbk4fdAcDPHCojPj4+io6OVlJSkm655ZbK5UlJSbr55pvPWz8wMFA7d+48Z9ns2bP17bffavny5YqMjKxmbMD5WctsenzZdn25I0sSD7sDgItxeNKzyZMna/To0YqJiVH//v317rvvKi0tTQkJCZIqTrFkZGRo0aJF8vDwUI8ePc7ZvkWLFvLz8ztvOeBKfjpTqvsWb9Hmoz/J29OiWbf30i19w8yOBQBOyeEyMmLECOXl5WnGjBnKyspSjx49tHLlSkVEREiSsrKyfnPOEcCVpZ8q1tgFm3T45Bk18vPSO6OjNaB9M7NjAYDTcnieETMwzwjqix3HT+uehZuVW1SqVkF+WjA+Vp1bNjI7FgCYoqrf3zybBqgh3+49oQeXpOhsmU1dQwO1cHw/hQRyxwwA/BbKCFADPtyYpmc+2ym7IQ3q2EyzR0WpkZ+32bEAoF6gjACXwDAMvfLNfr35XcXD7m6PDtPMW3ty6y4AOIAyAlRTabldU1bs0IqtGZKkR6/pqEnXduTWXQBwEGUEqIZCa5n+8MFWrT2YK08Pi/56Sw+N6NfG7FgAUC9RRgAHnSiwatyCzdqTVSB/H0+9NSpKV3VuYXYsAKi3KCOAAw6cKNS4BZuVcfqsmjX01YJx/dQzLMjsWABQr1FGgCradOSUJr6/WQXWcrVrHqD3x8cqPNjf7FgAUO9RRoAqWLkzS5MSt6m03K6oNo01b2w/NQnwMTsWALgEygjwGxasO6IZX6TKMKT4biF6/a6+8vP2NDsWALgMyghwEXa7oZdW7dU7PxyWJI2+LELPDe0uTw9u3QWAmkQZAS6gtNyuJ5dv1z+3ZUqSnrqhs/5wRXvmEAGAWkAZAX6h0FqmhA+Ste5gnrw8LJp1ey/dGhVmdiwAcFmUEeB//O8cIgE+nppzd7Qu79Tc7FgA4NIoI8B/Hcwp0tj5myrnEFk4vp96tGYOEQCobZQRQFLysZ804f3NOl1cpshmFXOItGnKHCIAUBcoI3B7Sakn9PDSrbKW2dU7vLHmj41R04a+ZscCALdBGYFbW7opTdM+3Sm7IV3VubneGhUlfx/+WgBAXWKvC7dkGIZe+88B/ePfByRJw2PC9NdbesrL08PkZADgfigjcDs2u6FnPtulpZvSJEkPX91Bk6/rxBwiAGASygjcirXMpkeWpuib1BOyWKQZQ7trdP+2ZscCALdGGYHbOF1cqonvb9GWYz/Jx8tDr9/ZRzf0CDU7FgC4PcoI3ELm6bMaO3+TDuQUKdDPS3PH9lNsZLDZsQAAoozADew/Uaix8zcpK9+qloF+ev+eWHVu2cjsWACA/6KMwKUlHzulexZuUf7ZMrVvHqBFE+LUunEDs2MBAP4HZQQuKyn1hB76cKtKyu3q26ax5o/tpyYBPmbHAgD8AmUELilxc5qmrqiYzOyaLi305sgoNfDxNDsWAOACKCNwKYZhaPb3h/Ty1/skSXdEh2nmrUxmBgDOjDICl2G3G5rxRaoWrj8qSXrgyvZ68vrOTGYGAE6OMgKXUFJu0+SPt+vLHVmSpOlDumn8wEiTUwEAqoIygnqvqKRc9y/eonUH8+TtadErw/toaO9WZscCAFQRZQT1Wm5RicYv2KydGfkK8PHUO6Nj9LuOzcyOBQBwAGUE9Vb6qWKNnrdRR/OK1TTARwvG91OvsMZmxwIAOIgygnopNbNAYxds0snCEoU1aaDFE+IU2SzA7FgAgGqgjKDe2Xg4TxPf36LCknJ1adlI798Tq5BAP7NjAQCqiTKCeuXr3dl6eGmKSsvtim0brPfGxiiogbfZsQAAl4Aygnrjo01p+uOnFbOqXtctRG/c1Vd+3syqCgD1HWUETu+Xs6qOiAnXX27pwayqAOAiKCNwana7oT9/maoF645KYlZVAHBFlBE4rTKbXU8u267PtmVKkv70+26a8DtmVQUAV0MZgVMqLi3XA0u26vt9J+XlYdHLd/TSLX3DzI4FAKgFlBE4ndPFpbpn4WZtTTstP28Pzbk7Wld1bmF2LABALaGMwKlk51s1Zv5G7T9RpEA/Ly0Y30/REcFmxwIA1CLKCJzG4ZNFGj1vkzJOn1VIoK8W3ROnzi0bmR0LAFDLKCNwCjuP52vcgk3KO1OqyGYBWnRPrMKD/c2OBQCoA5QRmG79oVzdtyhZRSXl6tE6UAvHx6pZQ1+zYwEA6ghlBKZatStLjyzdplKbXf3bNdW7Y6LVyI/p3QHAnVBGYJrEzWmauqJievfru4fotTuZ3h0A3BFlBKZ4e/UhvfjVXklM7w4A7o4ygjplGIZmfrVX7/5wWJKUcEV7PX0D07sDgDujjKDOlNvsmrpip5YlH5ckTbuxq+69vJ3JqQAAZqOMoE5Yy2x6ZGmKvkk9IU8Pi168tafuiAk3OxYAwAlQRlDrCq1lunfRFv14+JR8vDz05l19Fd+9pdmxAABOgjKCWpVXVKJxCzZrZ0a+Gvp66b0xMerfvqnZsQAAToQyglqTcfqsRs/dqMO5Z9Q0wEfv3xOrHq2DzI4FAHAylBHUioM5hRo9b5Oy8q1q3biBFk+IVbvmDc2OBQBwQpQR1Lgdx09r7PxN+qm4TB1aNNTiCbEKDWpgdiwAgJOijKBGrT+Uq3vf36IzpTb1DgvSgvGxCg7wMTsWAMCJUUZQY77ena2Hl6aotNyuAe2b6t0xMWroy68YAODX8U2BGrFsS7qe/mQHz5kBADiMMoJLNnfNYb3w5R5J0h3RYZp5a0+eMwMAqDLKCKrNMAy9mrRfb3x7UJI08XeRmnZTV54zAwBwSLX++Tp79mxFRkbKz89P0dHRWrNmzUXXXbFiha677jo1b95cgYGB6t+/v77++utqB4ZzsNsNTf98d2URefL6zhQRAEC1OFxGEhMTNWnSJE2bNk0pKSkaNGiQBg8erLS0tAuu/8MPP+i6667TypUrlZycrKuuukpDhgxRSkrKJYeHOcpsdj328TYt2nBMFov052E99OBVHSgiAIBqsRiGYTiyQVxcnKKiojRnzpzKZV27dtWwYcM0c+bMKv2M7t27a8SIEXr22WertH5BQYGCgoKUn5+vwMBAR+KihlnLbHpwyVb9Z2+OvDwsenVEHw3t3crsWAAAJ1TV72+HrhkpLS1VcnKypkyZcs7y+Ph4rV+/vko/w263q7CwUMHBwRddp6SkRCUlJZX/XVBQ4EhM1JICa5kmvr9Fm46ckq+Xh96+O1pXdWlhdiwAQD3n0Gma3Nxc2Ww2hYSEnLM8JCRE2dnZVfoZr7zyis6cOaPhw4dfdJ2ZM2cqKCio8hUezqPmzZZXVKKR7/2oTUdOqZGvlxZPiKOIAABqRLUuYP3ltQGGYVTpeoGlS5fqueeeU2Jiolq0uPgX2dSpU5Wfn1/5Sk9Pr05M1JCM02d1xzsbtCujQE0DfLT0vssUG3nxI1sAADjCodM0zZo1k6en53lHQXJycs47WvJLiYmJmjBhgpYtW6Zrr732V9f19fWVr6+vI9FQSw6dLNLouRuVyQPvAAC1xKEjIz4+PoqOjlZSUtI5y5OSkjRgwICLbrd06VKNGzdOH374oW666abqJUWd25WRr+Fvb1BmvlXtmwdoWUJ/iggAoMY5POnZ5MmTNXr0aMXExKh///569913lZaWpoSEBEkVp1gyMjK0aNEiSRVFZMyYMXrttdd02WWXVR5VadCggYKCgmrwo6AmbTpyShMWblZhSbl6tA7U++Nj1bQhR6sAADXP4TIyYsQI5eXlacaMGcrKylKPHj20cuVKRURESJKysrLOmXPknXfeUXl5uR588EE9+OCDlcvHjh2rhQsXXvonQI37bl+O/vBBsqxldsVGBmvu2BgF+nmbHQsA4KIcnmfEDMwzUnf+tT1TjyVuU7nd0NVdWmj2qCgeeAcAqJZamWcErm3ppjT98dOdMgxpaO9WemV4b3nzwDsAQC2jjECS9M7qQ5r51V5J0qi4Nppxcw95ejC9OwCg9lFG3JxhGHr5632a/f0hSdIfrmyvp67vzHNmAAB1hjLixux2Q89+vksf/FhxwfFTN3TWA1d2MDkVAMDdUEbcVJnNrieXbddn2zIrnrx7cw/dfVmE2bEAAG6IMuKGrGU2PfThVv17T8WTd18Z3ls392ltdiwAgJuijLiZopJy3fv+Fm04nCdfLw/NHhWla7r++lT+AADUJsqIGzldXKqxCzZre/ppBfh4au7YfurfvqnZsQAAbo4y4iZyCqwaPW+T9p0oVGN/b70/Pla9wxubHQsAAMqIO0g/VaxRczcq7VSxWjTy1QcT49QppJHZsQAAkEQZcXkHcwp199xNyi6wqk2wvz6YEKc2Tf3NjgUAQCXKiAvblZGvMfM36dSZUnVs0VAfTIxTSKCf2bEAADgHZcRFbTpyShMWblZhSbl6hQVp4fhYBQf4mB0LAIDzUEZc0Pf7cpTwQbKsZXbFRQZr7tgYNfLzNjsWAAAXRBlxMSt3ZunRj1JUZjN0VefmmnN3tPy8Pc2OBQDARVFGXMjHW9I15ZMdshvS73uF6tXhfeTj5WF2LAAAfhVlxEXMX3tEM75IlSTd2S9cf7mlpzw9ePIuAMD5UUbqOcMw9Ma3B/Vq0n5J0r2DIvXHG7vKYqGIAADqB8pIPWYYhv66co/eW3NEkjT5uk56+OoOFBEAQL1CGamnbHZDz3y2U0s3pUuSnv19N93zu0iTUwEA4DjKSD1UZrPrscRt+mJHljws0ou39dLwmHCzYwEAUC2UkXrGWmbTA0u26tu9OfL2tOgfI/rqpl6hZscCAKDaKCP1SFFJuSa+v1k/Hj4lP28PvX13tK7s3MLsWAAAXBLKSD1xurhUYxds1vb002ro66V5Y2MU166p2bEAALhklJF6IKfQqtFzN2nfiUI19vfWonti1SussdmxAACoEZQRJ3f8p2LdPXejjuYVq3kjX30wIU6dWzYyOxYAADWGMuLEDp8s0t1zNyoz36rWjRtoycQ4tW0WYHYsAABqFGXESe3JKtDoeRuVW1Sqds0DtGRinEKDGpgdCwCAGkcZcUIpaT9p7PxNKrCWq1tooBZNiFWzhr5mxwIAoFZQRpzM+kO5mvj+FhWX2hTVprEWjI9VUANvs2MBAFBrKCNO5Nu9J/SHD7aqpNyugR2a6t3RMQrw5X8RAMC18U3nJP61PVOPJW5Tud3QtV1D9ObIvvLz9jQ7FgAAtY4y4gQSN6dpyoqdMgzp5j6t9Lc7esvb08PsWAAA1AnKiMnmrT2iP3+RKkm6K7aNXhjWQ54eFpNTAQBQdygjJjEMQ298e1CvJu2XJN07KFJ/vLGrLBaKCADAvVBGTGAYhl78aq/e+eGwJOmxazvpkWs6UEQAAG6JMlLH7HZDf/rnLi3ZmCZJ+tPvu2nC7yJNTgUAgHkoI3Wo3GbXk8t36NOUDFks0sxbeurO2DZmxwIAwFSUkTpSUm7TI0tT9PXuE/LysOjVEX00tHcrs2MBAGA6ykgdKC4t1/2Lk7XmQK58vDw0e2SUru0WYnYsAACcAmWklhVYyzRh4WZtPvqT/H089d6YGA3s0MzsWAAAOA3KSC06daZUY+dv0s6MfDXy89LC8bGKjmhidiwAAJwKZaSW5BRYNWruRh3IKVLTAB8tmhCr7q2CzI4FAIDToYzUgvRTxbp73kYdyytWy0A/fTAxTh1aNDQ7FgAATokyUsMOnSzS3XM3KivfqvDgBvpw4mUKD/Y3OxYAAE6LMlKD9mQVaPS8jcotKlWHFg31wYQ4tQzyMzsWAABOjTJSQ1LSftLY+ZtUYC1X91aBWnRPrJo29DU7FgAATo8yUgN+PJynCQs360ypTVFtGmvB+FgFNfA2OxYAAPUCZeQSfbcvRwmLk1VSbteA9k313pgYBfgyrAAAVBXfmpdg5c4sPfpRispshq7t2kJvjoySn7en2bEAAKhXKCPVtDz5uJ5avl12Q/p9r1D9fUQfeXt6mB0LAIB6hzJSDYs3HNWf/rlbkjQiJlx/vbWnPD0sJqcCAKB+oow4aM73h/TSqr2SpPED2+pPN3WTB0UEAIBqo4xUkWEYeuWb/Xrzu4OSpIev7qDJ13WSxUIRAQDgUlBGqsBuNzTji1QtXH9UkjRlcBclXNHe3FAAALgIyshvsNkNTV2xQx9vOS5J+vPN3TW6f1tzQwEA4EIoI7+izGbXY4nb9MWOLHlYpJdv763bosPMjgUAgEuhjFyEtcymhz7cqn/vyZG3p0Wv39lXg3uGmh0LAACXQxm5gDMl5bpv8RatO5gnXy8PvT06Wld1bmF2LAAAXBJl5Bfyz5Zp/IJN2pp2WgE+npo3rp8ua9fU7FgAALgsysj/yCsq0Zj5m7Q7s0CBfl56/55Y9W3TxOxYAAC4tGrNXz579mxFRkbKz89P0dHRWrNmza+uv3r1akVHR8vPz0/t2rXT22+/Xa2wtSk736oR7/6o3ZkFatbQRx/d158iAgBAHXC4jCQmJmrSpEmaNm2aUlJSNGjQIA0ePFhpaWkXXP/IkSO68cYbNWjQIKWkpOiPf/yjHnnkEX3yySeXHL6mpJ8q1vB3NuhgTpFCg/yUeH9/dWsVaHYsAADcgsUwDMORDeLi4hQVFaU5c+ZULuvatauGDRummTNnnrf+008/rc8//1x79uypXJaQkKDt27drw4YNVXrPgoICBQUFKT8/X4GBNVsSDp0s0qj3Niq7wKo2wf5aMjFO4cH+NfoeAAC4o6p+fzt0ZKS0tFTJycmKj48/Z3l8fLzWr19/wW02bNhw3vrXX3+9tmzZorKysgtuU1JSooKCgnNetSE1s0Aj3tmg7AKrOrZoqGUJ/SkiAADUMYfKSG5urmw2m0JCQs5ZHhISouzs7Atuk52dfcH1y8vLlZube8FtZs6cqaCgoMpXeHi4IzGrxG439FjiNuUWlapH60Al3t9fIYF+Nf4+AADg11XrAtZfPhzOMIxffWDchda/0PKfTZ06Vfn5+ZWv9PT06sT8VR4eFr01qq+u6xaiD++9TMEBPjX+HgAA4Lc5dGtvs2bN5Onped5RkJycnPOOfvysZcuWF1zfy8tLTZteeP4OX19f+fr6OhKtWjq0aKT3xsTU+vsAAICLc+jIiI+Pj6Kjo5WUlHTO8qSkJA0YMOCC2/Tv3/+89b/55hvFxMTI29vbwbgAAMDVOHyaZvLkyZo7d67mz5+vPXv26LHHHlNaWpoSEhIkVZxiGTNmTOX6CQkJOnbsmCZPnqw9e/Zo/vz5mjdvnp544oma+xQAAKDecngG1hEjRigvL08zZsxQVlaWevTooZUrVyoiIkKSlJWVdc6cI5GRkVq5cqUee+wxvfXWW2rVqpVef/113XbbbTX3KQAAQL3l8DwjZqjNeUYAAEDtqJV5RgAAAGoaZQQAAJiKMgIAAExFGQEAAKaijAAAAFNRRgAAgKkoIwAAwFSUEQAAYCrKCAAAMJXD08Gb4edJYgsKCkxOAgAAqurn7+3fmuy9XpSRwsJCSVJ4eLjJSQAAgKMKCwsVFBR00T+vF8+msdvtyszMVKNGjWSxWGrs5xYUFCg8PFzp6ek886aKGDPHMF6OY8wcx5g5hvFyXHXHzDAMFRYWqlWrVvLwuPiVIfXiyIiHh4fCwsJq7ecHBgbyC+kgxswxjJfjGDPHMWaOYbwcV50x+7UjIj/jAlYAAGAqyggAADCVW5cRX19fTZ8+Xb6+vmZHqTcYM8cwXo5jzBzHmDmG8XJcbY9ZvbiAFQAAuC63PjICAADMRxkBAACmoowAAABTUUYAAICpXL6MzJ49W5GRkfLz81N0dLTWrFnzq+uvXr1a0dHR8vPzU7t27fT222/XUVLn4ciYrVixQtddd52aN2+uwMBA9e/fX19//XUdpjWfo79jP1u3bp28vLzUp0+f2g3ohBwds5KSEk2bNk0RERHy9fVV+/btNX/+/DpKaz5Hx2vJkiXq3bu3/P39FRoaqvHjxysvL6+O0prvhx9+0JAhQ9SqVStZLBZ99tlnv7mNO+/7HR2vWtnvGy7so48+Mry9vY333nvPSE1NNR599FEjICDAOHbs2AXXP3z4sOHv7288+uijRmpqqvHee+8Z3t7exvLly+s4uXkcHbNHH33UeOmll4xNmzYZ+/fvN6ZOnWp4e3sbW7durePk5nB0vH52+vRpo127dkZ8fLzRu3fvugnrJKozZkOHDjXi4uKMpKQk48iRI8bGjRuNdevW1WFq8zg6XmvWrDE8PDyM1157zTh8+LCxZs0ao3v37sawYcPqOLl5Vq5caUybNs345JNPDEnGp59++qvru/u+39Hxqo39vkuXkdjYWCMhIeGcZV26dDGmTJlywfWfeuopo0uXLucsu//++43LLrus1jI6G0fH7EK6detmPP/88zUdzSlVd7xGjBhhPPPMM8b06dPdrow4OmZfffWVERQUZOTl5dVFPKfj6Hi9/PLLRrt27c5Z9vrrrxthYWG1ltGZVeXLlX3//6vKeF3Ipe73XfY0TWlpqZKTkxUfH3/O8vj4eK1fv/6C22zYsOG89a+//npt2bJFZWVltZbVWVRnzH7JbrersLBQwcHBtRHRqVR3vBYsWKBDhw5p+vTptR3R6VRnzD7//HPFxMRo1qxZat26tTp16qQnnnhCZ8+erYvIpqrOeA0YMEDHjx/XypUrZRiGTpw4oeXLl+umm26qi8j1krvv+y9VTez368WD8qojNzdXNptNISEh5ywPCQlRdnb2BbfJzs6+4Prl5eXKzc1VaGhoreV1BtUZs1965ZVXdObMGQ0fPrw2IjqV6ozXgQMHNGXKFK1Zs0ZeXi771++iqjNmhw8f1tq1a+Xn56dPP/1Uubm5euCBB3Tq1CmXv26kOuM1YMAALVmyRCNGjJDValV5ebmGDh2qN954oy4i10vuvu+/VDWx33fZIyM/s1gs5/y3YRjnLfut9S+03JU5OmY/W7p0qZ577jklJiaqRYsWtRXP6VR1vGw2m0aOHKnnn39enTp1qqt4TsmR3zG73S6LxaIlS5YoNjZWN954o1599VUtXLjQLY6OSI6NV2pqqh555BE9++yzSk5O1qpVq3TkyBElJCTURdR6i31/9dTUft9l/2nWrFkzeXp6nvevh5ycnPMa8M9atmx5wfW9vLzUtGnTWsvqLKozZj9LTEzUhAkTtGzZMl177bW1GdNpODpehYWF2rJli1JSUvTQQw9JqviiNQxDXl5e+uabb3T11VfXSXazVOd3LDQ0VK1btz7nMeRdu3aVYRg6fvy4OnbsWKuZzVSd8Zo5c6YGDhyoJ598UpLUq1cvBQQEaNCgQXrhhRf4V/4FuPu+v7pqcr/vskdGfHx8FB0draSkpHOWJyUlacCAARfcpn///uet/8033ygmJkbe3t61ltVZVGfMpIpmPG7cOH344YdudV7a0fEKDAzUzp07tW3btspXQkKCOnfurG3btikuLq6uopumOr9jAwcOVGZmpoqKiiqX7d+/Xx4eHgoLC6vVvGarzngVFxfLw+PcXbunp6ek///XPs7l7vv+6qjx/X61L32tB36+JW7evHlGamqqMWnSJCMgIMA4evSoYRiGMWXKFGP06NGV6/98e9djjz1mpKamGvPmzXOr27sMw/Ex+/DDDw0vLy/jrbfeMrKysipfp0+fNusj1ClHx+uX3PFuGkfHrLCw0AgLCzNuv/12Y/fu3cbq1auNjh07GhMnTjTrI9QpR8drwYIFhpeXlzF79mzj0KFDxtq1a42YmBgjNjbWrI9Q5woLC42UlBQjJSXFkGS8+uqrRkpKSuXt0Oz7z+XoeNXGft+ly4hhGMZbb71lREREGD4+PkZUVJSxevXqyj8bO3asccUVV5yz/vfff2/07dvX8PHxMdq2bWvMmTOnjhObz5Exu+KKKwxJ573Gjh1b98FN4ujv2P9yxzJiGI6P2Z49e4xrr73WaNCggREWFmZMnjzZKC4uruPU5nF0vF5//XWjW7duRoMGDYzQ0FBj1KhRxvHjx+s4tXm+++67X90vse8/l6PjVRv7fYthcNwOAACYx2WvGQEAAPUDZQQAAJiKMgIAAExFGQEAAKaijAAAAFNRRgAAgKkoIwAAwFSUEQAAYCrKCAAAMBVlBAAAmIoyAgAATEUZAQAApvo/zLmamWLFPV4AAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "a = start[0]\n", + "b = end[0]\n", + "domain = np.linspace(a,b, 100)\n", + "plt.plot(domain, np.sin(domain))\n", + "#plt.show()" + ] + }, + { + "cell_type": "code", + "execution_count": 5, + "id": "d8fcfba3", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Integral in first domain: 0.6173165676349102\n" + ] + } + ], + "source": [ + "print('Integral in first domain: {}'.format(sin_integral(a,b)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "47ed9e21", + "metadata": {}, + "outputs": [], + "source": [ + "a = start[1]\n", + "b = end[1]\n", + "domain = np.linspace(a,b, 100)\n", + "plt.plot(domain, np.sin(domain))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf914ead", + "metadata": {}, + "outputs": [], + "source": [ + "print('Integral in second domain: {}'.format(sin_integral(a,b)))" + ] + }, + { + "cell_type": "markdown", + "id": "387e2fd6", + "metadata": {}, + "source": [ + "#### $2^{nd}$ Domain: $[\\pi, \\frac{5\\pi}{4}]$\n", + "\n", + "In this domain, the function is not strictly positively defined and the integral is negative!" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "9bc7ab74", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "[]" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAi8AAAGdCAYAAADaPpOnAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABCUklEQVR4nO3deVxU5eIG8OfMDDAgMCqbC4i74q6o7HorxdyxTE0FtVxI0dTqXs1umXUv93ZvXRNFza1yQS21XMgkK2URVMQV3EAFFURQZxBknfP7Q+UXiQrKzJnl+X4+84eHc4aHc7nM0znveV9BFEURREREREZCJnUAIiIiotpgeSEiIiKjwvJCRERERoXlhYiIiIwKywsREREZFZYXIiIiMiosL0RERGRUWF6IiIjIqCikDlDXtFotrl+/Djs7OwiCIHUcIiIiqgFRFFFQUIAmTZpAJnvytRWTKy/Xr1+Hm5ub1DGIiIjoGWRlZcHV1fWJ+5hcebGzswNw/4e3t7eXOA0RERHVhEajgZubW+Xn+JOYXHl5eKvI3t6e5YWIiMjI1GTIBwfsEhERkVFheSEiIiKjwvJCRERERoXlhYiIiIwKywsREREZFZYXIiIiMiosL0RERGRUWF6IiIjIqLC8EBERkVFheSEiIiKjwvJCRERERoXlhYiIiIwKy0stLNx5Bl8dTEdpuVbqKERERGbL5FaV1pXU6xp8nXAZALApKRMLBndAPw/nGq1+SURERHWHV15qqH0jO3z2ahc42lrhcn4Rpnx7FMFrDuNcToHU0YiIiMyKIIqiKHWIuqTRaKBSqaBWq2Fvb1/n73+3pBzLfruINbGXUFqhhUwAxnm5Y27/tmhQz7LOvx8REZE5qM3nN8vLM8rML0L4T2n46XQOAEBlbYG5/dtinFczKOS8oEVERFQbLC96KC8PJaTnYdGuVJx9cPuorYstPhraEX6tHXX+vYmIiEwFy4seywsAlFdoEXUkC5/vO4c7RWUAgJc7NsKCwR5wa2ijlwxERETGjOVFz+XloTtFpVj8ywWsT7yCCq0IK4UM0/q2wlt9W8HaUq7XLERERMaE5UWi8vLQ2RwNPt6ZikMZ+QCAJiolFgzugEGdG/HRaiIiomqwvEhcXgBAFEX8dDoH/9iThmt37gEA/Fo7YOHQjmjjYidZLiIiIkPE8mIA5eWhe6UVWH4gHSsO3J+ZVyETMNG3Od7u1wZ2Sgup4xERERmE2nx+85leHbO2lGNu/7b4ZU5f9O/ggnKtiNVxl/DCfw9gR8pVmFh3JCIi0jleedGz38/lYtGuVGTkFQIAejdviI+Hd4RHY8PLSkREpC+8bWTA5QUASsorsCbuEiL2X8S9sgrIZQKCvd0xN7At7HkriYiIzBBvGxk4K4Uc0//SGr+80xeDOjdChVbE1wmX8SJvJRERET0Vr7wYgNgLN/HRj2cqbyV5tWiIT4I6oS2fSiIiIjPBKy9GJqCNE36aHYD3BrSD0kKGpEu3MOjLWPwzOg2FJeVSxyMiIjIoLC8Gwkohx4wXWiPmD08lfXUwA/2+OIC9p7N5K4mIiOgBlhcD49bQBqtCemLtxJ5wa2iNbHUxQjccwxtfH0FmfpHU8YiIiCTH8mKgXmzvgn2z+yLshdawkAv47dxN9P/fAUTsv4CS8gqp4xEREUmG5cWAWVvK8e6Adtg7uw98WzmgpFyLz2POY9CXsTiUni91PCIiIkmwvBiBVk622DjZC1+O6QZHW0uk3yzE66sSMXfrceTfLZE6HhERkV6xvBgJQRAwvFtT7J/7F4z3bgZBALYfu4YXPz+ALUcyodVyQC8REZkHlhcjo7KxwKdBnbH9LV90aGwP9b0y/G3bKYz5KhEXbhRIHY+IiEjnWF6MVPdmDbAzzA8fDPaAtYUchy/fwqAlsfjvz+dQXMYBvUREZLpYXoyYQi7D5ICW+OWdvujn4YyyChFLf7uIAYsPIu5CntTxiIiIdILlxQQ0rW+NVSE9sWK8JxrZK3Elvwjj1yRh7hYO6CUiItPD8mIiBEHAy50aIWZuH0zwcb8/oDflGvp9cQDfJ3OxRyIiMh0sLybGTmmBj4d3wva3fNG+kR1uF5Xh3e9OYPyaJFx+sPAjERGRMWN5MVHdmzXArpn+mDewPZQWMsRfzMeAxQex/Pd0lFVopY5HRET0zPRSXiIjI9GiRQsolUp4enoiNjb2ifuXlJRgwYIFcHd3h5WVFVq1aoW1a9fqI6pJsZDLENq3FX6e3Qf+rR1RUq7Fv/eexbCl8TiRdUfqeERERM9E5+Vly5YtmD17NhYsWICUlBQEBARg4MCByMzMfOwxo0aNwv79+7FmzRqcO3cOUVFRaN++va6jmix3h3pY/2ZvfP5aVzSwsUBatgYjIuPxye5UFJWWSx2PiIioVgRRxyM5vby80KNHDyxfvrxym4eHB4KCghAeHv7I/nv37sWYMWOQkZGBhg0b1vr7aTQaqFQqqNVq2NvbP1d2U5R/twSf7E7FD8evAwBcG1jjnyM6o09bJ4mTERGROavN57dOr7yUlpYiOTkZgYGBVbYHBgYiISGh2mN27tyJnj174rPPPkPTpk3Rtm1bvPvuu7h3754uo5oNB1srLB7THesm9ULT+ta4evseQtYextytx3G7sFTqeERERE+l0OWb5+XloaKiAi4uLlW2u7i4ICcnp9pjMjIyEBcXB6VSiR07diAvLw/Tp0/HrVu3qh33UlJSgpKS/5/LRKPR1O0PYaJeaOeMfXP64D8/n8M3hy5j+7FrOHDuJj4e3hGDOzeGIAhSRyQiIqqWXgbs/vmDUBTFx344arVaCIKAjRs3onfv3hg0aBC++OILfP3119VefQkPD4dKpap8ubm56eRnMEX1rBRYOKwjtr3li7YutsgvLEXYphRMXZ+MG5piqeMRERFVS6flxdHREXK5/JGrLLm5uY9cjXmocePGaNq0KVQqVeU2Dw8PiKKIq1evPrL//PnzoVarK19ZWVl1+0OYgR4PHqt++6U2sJALiEm9gX5fHEDU4UxObkdERAZHp+XF0tISnp6eiImJqbI9JiYGvr6+1R7j5+eH69ev4+7du5Xbzp8/D5lMBldX10f2t7Kygr29fZUX1Z6VQo45/dti10x/dHVVoaC4HPO3n8K41UnIzC+SOh4REVElnd82mjt3LlavXo21a9ciLS0Nc+bMQWZmJkJDQwHcv3ISEhJSuf/YsWPh4OCASZMmITU1FQcPHsR7772HN954A9bW1rqOa/baN7LH9un3V6tWWsiQkH5/crs1cZdQoeVVGCIikp5OB+wCwOjRo5Gfn49FixYhOzsbnTp1QnR0NNzd3QEA2dnZVeZ8sbW1RUxMDGbOnImePXvCwcEBo0aNwqeffqrrqPSAXCZgckBL9PNwwbztJ5GYcQuf7E7FnpPX8dnILmjtbCd1RCIiMmM6n+dF3zjPS93SakVEHclEePRZ3C0ph6VChtn92mBqQEso5FxdgoiI6obBzPNCxk8mEzDOyx375vTBX9o5obRci8/2nsOIyASczeFj6UREpH8sL1QjTepbY93EXvjva11hr1Tg1DU1hkbE4ctfLnChRyIi0iuWF6oxQRAw0tMVMXP7op+HC8oqRPzvl/MYvjQeZ66rpY5HRERmguWFas3FXolVIZ74ckw3NLCxQGq2BsOXxuN/MedRWs6rMEREpFssL/RMBEHA8G5NsW9OX7zcsRHKtSK+3H8Bw5bG4fQ1XoUhIiLdYXmh5+JkZ4Xl43sg4vXuaGBjgbM5BQhaFo8veBWGiIh0hOWFnpsgCBjatQli5vbFwE73r8Is2X8Bw5fFI/U6n0giIqK6xfJCdcbR1gqR4/7/KkxatgbDlvKJJCIiqlssL1SnHl6F+eNYmP/9ch4jIuNxLqdA6nhERGQCWF5IJx6OhVnyenfUt7HA6WsaDI2IQ+TvF1HOqzBERPQcWF5IZwRBwLCuTbBvTh/083BGacX92XlHrjiEi7l3n/4GRERE1WB5IZ1ztlNiVUhP/Pe1rrBTKnA86w4GL4nFmrhL0HKlaiIiqiWWF9KLh7Pz7pvTB33aOqGkXItPdqfi9VWJyLpVJHU8IiIyIiwvpFeNVdb4ZlIv/HNEZ9hYypF06RYGLD6ITUmZMLEFzomISEdYXkjvBEHAWK9m2Pt2H/Ru3hBFpRV4f8cpTFx3BDnqYqnjERGRgWN5Ick0c7DB5qne+GCwBywVMhw4fxMDFh/EzhPXpY5GREQGjOWFJCWTCZgc0BJ7Zvqjc1MV1PfKMCsqBTM2HcPtwlKp4xERkQFieSGD0MbFDtun+2J2vzaQywTsOZmNwMUH8dvZXKmjERGRgWF5IYNhIZdhdr+22DHdF62dbXGzoASTvj6C+dtPobCkXOp4RERkIFheyOB0ca2P3TP98YZfCwBA1OFMDFoSi+QrtyRORkREhoDlhQyS0kKOD4d2wKbJXmiiUuJKfhFeW3EIn+09i9JyLi9ARGTOWF7IoPm2dsRPs/vglR5NoRWByN/TEbQsHudvcJFHIiJzxfJCBk9lbYEvRnXD8nE90MDGAqnZGgyJiMPq2AwuL0BEZIZYXshoDOzcGD/P7oO/tHNCabkWn+5Jw/g1Sbh+557U0YiISI9YXsioONsrsW5iL3wa1AnWFnIkpOfjZU5sR0RkVlheyOgIgoDx3u7YM8sfXd3qQ1NcjllRKXh7cwrU98qkjkdERDrG8kJGq6WTLb4P9cHbL92f2O7H49cxcPFBHErPlzoaERHpEMsLGTULuQxz+rfF96E+aO5gg+vqYoxdnYjw6DSUlFdIHY+IiHSA5YVMQvdmDbBnVgBe7+0GUQRWHszAiGUJuMBHqomITA7LC5mMelYKhL/SBV8Fe6JhPcvKR6q/jr8EUeQj1UREpoLlhUxOYMdG2Ds7AH3bOqGkXIuFu1Ixcd0R5BYUSx2NiIjqAMsLmSRnOyW+ntQLC4d2gJVChgPnb+LlxbH4JfWG1NGIiOg5sbyQyRIEARP9WmDXTH+0b2SHW4WlmPztUSzYcQr3SjmYl4jIWLG8kMlr62KHH8P8MCXg/irVG5MyMTgiFqevqSVORkREz4LlhcyClUKOBYM7YMObXnCxt0LGzUKMiIzHVwfTuT4SEZGRYXkhs+LfxhF73+6DAR1dUFYh4p/RZxGy9jBuaDiYl4jIWLC8kNlpUM8SK8Z7IvyVzrC2kCPuYh4GLD6In8/kSB2NiIhqgOWFzJIgCHi9dzPsnuWPTk3tcaeoDNPWJ3MwLxGREWB5IbPWyskW29/yw7Q+LQHcH8w7dGkczlznYF4iIkPF8kJmz1Ihw/xBHtjwphec7axwMfcuRixLwJo4zsxLRGSIWF6IHvBv44if3g5APw9nlFZo8cnuVEz6+ghuFpRIHY2IiP6A5YXoDxxsrbAqpCc+Gd4RVgoZfj93EwO/jMWB8zeljkZERA+wvBD9iSAICPZpjp1h/mjrYou8uyWYsPYwPt2dipJyDuYlIpIaywvRY7RrZIedYf4I8XEHAKyOu4RXlycg4+ZdiZMREZk3lheiJ1BayLFoeCesCumJBjYWOH1NgyERcdiWfJWDeYmIJMLyQlQD/Tu44Ke3+8C7ZUMUlVbgne9OYM6W4ygoLpM6GhGR2WF5IaqhRiolNk72xruBbSGXCfjh+HUMXhKHE1l3pI5GRGRWWF6IakEuExD2YhtsneaNpvWtkXmrCK8uT8DKA1zgkYhIX1heiJ6Bp3tDRL8dgMGdG6NcKyL8p7OYyDlhiIj0guWF6BmprC2wdGx3hL/SGUoLGQ6ev4mBXx7EQc4JQ0SkUywvRM/h4QKPu8L80c7FDnl3SxGy9jD+vfcsyiq0UscjIjJJeikvkZGRaNGiBZRKJTw9PREbG1uj4+Lj46FQKNCtWzfdBiR6Tm1c7PBjmB/GezcDACz/PR2jVh5C1q0iiZMREZkenZeXLVu2YPbs2ViwYAFSUlIQEBCAgQMHIjMz84nHqdVqhISE4KWXXtJ1RKI6obSQ49Ogzlg+rgfslAqkZN7BoCWxiD6VLXU0IiKTIog6nmnLy8sLPXr0wPLlyyu3eXh4ICgoCOHh4Y89bsyYMWjTpg3kcjl++OEHHD9+vEbfT6PRQKVSQa1Ww97e/nnjEz2TrFtFmLU5BSmZdwAA47ya4e9DOkBpIZc2GBGRgarN57dOr7yUlpYiOTkZgYGBVbYHBgYiISHhscetW7cO6enp+Oijj576PUpKSqDRaKq8iKTm1tAGW6f5ILRvKwDAxqRMBC2Lx8VcLi1ARPS8dFpe8vLyUFFRARcXlyrbXVxckJOTU+0xFy5cwLx587Bx40YoFIqnfo/w8HCoVKrKl5ubW51kJ3peFnIZ5g1sj2/e6A2HepY4m1OAoRFx+D75qtTRiIiMml4G7AqCUOXfoig+sg0AKioqMHbsWHz88cdo27Ztjd57/vz5UKvVla+srKw6yUxUV/q2dcJPbwfAr7UD7pVV4N3vTmDuluMoLCmXOhoRkVF6+qWN5+Do6Ai5XP7IVZbc3NxHrsYAQEFBAY4ePYqUlBSEhYUBALRaLURRhEKhwL59+/Diiy9WOcbKygpWVla6+yGI6oCzvRLfvuGF5b9fxBcx57E95RqOX72Dpa/3QIcmHJtFRFQbOr3yYmlpCU9PT8TExFTZHhMTA19f30f2t7e3x6lTp3D8+PHKV2hoKNq1a4fjx4/Dy8tLl3GJdOrh0gKbp/qgkb0SGTcLERQZj41JV7hCNRFRLej0ygsAzJ07F8HBwejZsyd8fHzw1VdfITMzE6GhoQDu3/a5du0avv32W8hkMnTq1KnK8c7OzlAqlY9sJzJWvVvcX1rg3e9O4NezuViw4zQS0vMR/kpn2CstpI5HRGTwdF5eRo8ejfz8fCxatAjZ2dno1KkToqOj4e7uDgDIzs5+6pwvRKamYT1LrA7piTVxl/DvvWex52Q2Tl1VY9nYHujsqpI6HhGRQdP5PC/6xnleyNikZN5G2KYUXLtzD5ZyGd4f1B4TfJtXO6idiMhUGcw8L0T0dN2bNUD0rAAEdnBBaYUWC3elInRDMtRFZVJHIyIySCwvRAZAZWOBlcGe+GhoB1jIBfx85gYGR8TieNYdqaMRERkclhciAyEIAib5tcC2t3zh1tAaV2/fw2srErA27hKfRiIi+gOWFyID08W1PnbPDMDATo1QViFi0e5UTFvP20hERA+xvBAZIJW1BSLH9cDHwzrCUi7DvlTeRiIieojlhchACYKACb7Nse0tXzRraMPbSERED7C8EBm4zq4q7J7lX+U20lsbjkF9j7eRiMg8sbwQGQF75f3bSAsfPI2090wOhkTE4uTVO1JHIyLSO5YXIiMhCAIm+rXAd6G+cG1gjaxb9zBy+SF8e+gybyMRkVlheSEyMt3c6mPPzAD0fzCp3Yc/nkFYVAoKinkbiYjMA8sLkRFS2Vjgq2BPfDDYAwqZgD0nszFsaTxSr2ukjkZEpHMsL0RGShAETA5oiS3TfNBEpcSlvEKMiIzH5sOZvI1ERCaN5YXIyHm6N8CeWQF4oZ0TSsq1mLf9FN7ZegJFpeVSRyMi0gmWFyIT0KCeJdZM6IW/vtwOMgHYnnINw5fG42JugdTRiIjqHMsLkYmQyQRM/0trbJriDWc7K1zIvYthS+PxQ8o1qaMREdUplhciE+Pd0gF7ZgXAt5UDikorMHvLcSzYcQrFZRVSRyMiqhMsL0QmyMnOCuvf9MKsF1tDEICNSZkYuSIBmflFUkcjInpuLC9EJkouEzA3sB3WTeyFBjYWOH1Ng8ERsYhJvSF1NCKi58LyQmTi/tLOGXtmBaB7s/ooKC7HlG+PIjw6DeUVWqmjERE9E5YXIjPQpL41tkz1wRt+LQAAKw9mYOyqJORqiiVORkRUeywvRGbCUiHDh0M7IHJcD9haKXD48i0MWhKHhPQ8qaMREdUKywuRmRnUuTF2hvmhfSM75N0twfjVSVj220VotZyVl4iMA8sLkRlq6WSLHdP9MNLTFVoR+M/P5zDl26NQF3FxRyIyfCwvRGbK2lKO/77WFf9+tTOsFDLsP5uLwRGxOHVVLXU0IqInYnkhMnOjezXD9um+aNbQBldv38OryxOwKYmLOxKR4WJ5ISJ0bKLCrpn+6OfhgtIKLd7fcQrvfHcC90o5Ky8RGR6WFyICAKisLbAqxBPzBraHXCZg+7FrGBEZj0t5hVJHIyKqguWFiCoJgoDQvq2w4U0vONpa4WxOAYZFxGHv6WypoxERVWJ5IaJH+LRyQPQsf/Ru3hAFJeUI3XAM/+SsvERkIFheiKhazvZKbJzihal9WgIAvjqYgbGrOSsvEUmP5YWIHstCLsP7gzyw/OGsvJduYXBEHA5fuiV1NCIyYywvRPRUAx/MytvOxQ43C0rw+qpErI7N4OPURCQJlhciqpGWTrbYMcMXQd2aoEIr4tM9aZi+8RgKijkrLxHpF8sLEdWYjaUC/xvdDZ8EdYKFXMBPp3MwfFk8zt8okDoaEZkRlhciqhVBEBDs7Y6t03zQWKVExs1CDF8aj50nrksdjYjMBMsLET2T7s0aYPdMf/i3dsS9sgrMikrBwp1nUFrOx6mJSLdYXojomTnYWuGbN3pjxgutAABfJ1zG66sScYOPUxORDrG8ENFzkcsEvDegPVaF9ISdUoHkK7cxeEkcEjPypY5GRCaK5YWI6kT/Di7YFeaP9o3skHe3BONWJ2HVQT5OTUR1j+WFiOpMc8d62DHdDyO6N0WFVsQ/otMwY9Mx3C0plzoaEZkQlhciqlPWlnJ8MaorPhneERZyAdGncjB8aRwu5vJxaiKqGywvRFTnBEFAsE9zbJnmg0b2SqQ/eJx6z0muTk1Ez4/lhYh0pkezBtg9yx8+LR1QWFqBGZuO4R97Urk6NRE9F5YXItIpR1srrH+zN6b1vb869arYSxi3Ogk3C0okTkZExorlhYh0TiGXYf5AD6wY3wP1LOVIunQLQyPicCzzttTRiMgIsbwQkd683KkxfgzzRyunesjRFGP0ykNYn3iFj1MTUa2wvBCRXrV2tsWPYf4Y1LkRyipE/P2H03j3u5MoLquQOhoRGQmWFyLSO1srBZaN7YH3B7WHTAC2HbuKVyITkHWrSOpoRGQEWF6ISBKCIGBqn1bYMNkLDvUskZqtwZCIOBw4f1PqaERk4FheiEhSvq0csWumP7q61Yf6XhkmrjuMiP0XoNVyHAwRVY/lhYgk16S+NbZO88ZYr2YQReDzmPOYuv4oNMVlUkcjIgOkl/ISGRmJFi1aQKlUwtPTE7GxsY/dd/v27ejfvz+cnJxgb28PHx8f/Pzzz/qISUQSslLI8c8RnfHZq11gqZDhl7RcDF8aj/M3uKwAEVWl8/KyZcsWzJ49GwsWLEBKSgoCAgIwcOBAZGZmVrv/wYMH0b9/f0RHRyM5ORkvvPAChg4dipSUFF1HJSIDMKqXG74P9UHT+ta4lFeIoGXx2H3yutSxiMiACKKOJ1jw8vJCjx49sHz58sptHh4eCAoKQnh4eI3eo2PHjhg9ejQ+/PDDp+6r0WigUqmgVqthb2//zLmJSFq3CksxM+oY4i/mAwCm9mmJvw5oB4Wcd7uJTFFtPr91+legtLQUycnJCAwMrLI9MDAQCQkJNXoPrVaLgoICNGzYsNqvl5SUQKPRVHkRkfFrWM8S30zqjdC+rQAAXx3MQMjaw8i/y2UFiMydTstLXl4eKioq4OLiUmW7i4sLcnJyavQen3/+OQoLCzFq1Khqvx4eHg6VSlX5cnNze+7cRGQYFHIZ5g1sj+Xj7i8rkJCej6ERcTh59Y7U0YhIQnq5/ioIQpV/i6L4yLbqREVFYeHChdiyZQucnZ2r3Wf+/PlQq9WVr6ysrDrJTESGY2Dnxvhhhh9aOtbDdXUxRq44hK1H+P91InOl0/Li6OgIuVz+yFWW3NzcR67G/NmWLVvw5ptvYuvWrejXr99j97OysoK9vX2VFxGZnjYudvghzA/9O7igtFyLv247ifd3nEJJOZcVIDI3Oi0vlpaW8PT0RExMTJXtMTEx8PX1fexxUVFRmDhxIjZt2oTBgwfrMiIRGRF7pQVWjvfEu4FtIQjApqRMjPkqETnqYqmjEZEe6fy20dy5c7F69WqsXbsWaWlpmDNnDjIzMxEaGgrg/m2fkJCQyv2joqIQEhKCzz//HN7e3sjJyUFOTg7UarWuoxKREZDJBIS92AZrJ/aCvVKBlMw7GBIRi6SMfKmjEZGe6Ly8jB49GosXL8aiRYvQrVs3HDx4ENHR0XB3dwcAZGdnV5nzZeXKlSgvL8eMGTPQuHHjytfbb7+t66hEZEReaOeMXTP90b6RHfLulmLc6iSsi78EHc/+QEQGQOfzvOgb53khMi9FpeWYt+0Udp64P5HdiO5N8c8RnWFtKZc4GRHVhsHM80JEpGs2lgp8OaYbPhjsAblMwI6Ua3h1eQKybhVJHY2IdITlhYiMniAImBzQEhve9IJDPUukZmswdGkcYi/clDoaEekAywsRmQyfVg7YNdMfXVxVuFNUhglrD2PFgXSOgyEyMSwvRGRSmtS3xtZpPnjN0xVaEfjXT2cRFpWCwpJyqaMRUR1heSEik6O0kOOzkV3wSVAnWMgF7DmZjVciE3A5r1DqaERUB1heiMgkCYKAYG93RE3xhpOdFc7dKMCwpXH47Wyu1NGI6DmxvBCRSevZvCF2z/RHj2b1oSkuxxvfHEHE/gvQajkOhshYsbwQkclzsVdi81QfjPNqBlEEPo85j9ANySgoLpM6GhE9A5YXIjILlgoZ/jGiM/79amdYymXYl3oDQcvikX7zrtTRiKiWWF6IyKyM7tUMW6Z5o5G9Euk3CxG0NB4xqTekjkVEtcDyQkRmp3uzBtg10x+9mzdEQUk5pnx7FP+LOc9xMERGguWFiMySk50VNk7xwkTf5gCAL/dfwNT1R6HhOBgig8fyQkRmy0Iuw8JhHfHf17rCUiHDL2m5CFoWj4u5HAdDZMhYXojI7I30dMX3oT5orFIi42YhgpbF4+czOVLHIqLHYHkhIgLQxbX+/XEwLRribkk5pq1Pxhf7znEcDJEBYnkhInrA0dYKGyf//ziYJb9exJRvOQ6GyNCwvBAR/cHDcTCfv9YVVgoZ9p/NRdDSeFy4USB1NCJ6gOWFiKgar3q64vtQXzRRKZGRx3EwRIaE5YWI6DE6u6qwa6Y/vFs2RGFpBcfBEBkIlhcioidwsLXChje98IZfCwD3x8FM/vYo1Pc4DoZIKiwvRERPoZDL8OHQDvhi1P1xML+efTgfDMfBEEmB5YWIqIZe6fH/42Au5RUiaFkCx8EQSYDlhYioFjq7qrBzpj+8/jgfDNdFItIrlhciolpytLXChj/OB8N1kYj0iuWFiOgZ/HE+mMp1kZZyXSQifWB5ISJ6Dq/+cV2kB/PBxKTekDoWkUljeSEiek5/XhdpyrdH8T+OgyHSGZYXIqI68HBdpAk+7gCAL/dfwNT1ySjgOBiiOsfyQkRURyzkMnw8vBM+G9kFlnIZfkm7gaBl8Ui/yXEwRHWJ5YWIqI6N6umGraE+aGSvRPrNQgQtjcevZzkOhqiusLwQEelAN7f62DnTD72aN0BBSTne/OYoIvZf4DgYojrA8kJEpCPOdkpsnOyN8d7NIIrA5zHnMX3jMdwtKZc6GpFRY3khItIhS4UMnwZ1xr9e6QwLuYC9Z3LwSmQ8LucVSh2NyGixvBAR6cGY3s2weaoPnO2scP7GXQxbGocD529KHYvIKLG8EBHpiad7A+ya6Y/uzepDU1yOSesOY/nv6RBFjoMhqg2WFyIiPXKxV2LzVG+M6eUGrQj8e+9ZhEWloKiU42CIaorlhYhIz6wUcoS/0hmfBnWCQiZgz8lsvBKZgKxbRVJHIzIKLC9ERBIQBAHjvd2xaYo3HG0tcTanAEOXxiH+Yp7U0YgMHssLEZGEerdoiF0z/dHFVYU7RWUIXpOE1bEZHAdD9AQsL0REEmusssbWaT54tYcrtCLw6Z40zN16AsVlFVJHIzJILC9ERAZAaSHHf1/rgo+GdoBcJmBHyjWMXJGAa3fuSR2NyOCwvBARGQhBEDDJrwXWv9kbDWwscPqaBsMi4pCUkS91NCKDwvJCRGRgfFs5YmeYPzo0tkd+YSnGrU7Ct4cucxwM0QMsL0REBsitoQ22veWLYV2boFwr4sMfz+Bv206ipJzjYIhYXoiIDJS1pRxfjumG+QPbQyYAW49exeiVibihKZY6GpGkWF6IiAyYIAiY1rcVvp7UGyprCxzPuoMhEXFIvnJb6mhEkmF5ISIyAn3aOmFnmB/audjhZkEJxnx1CFGHM6WORSQJlhciIiPh7lAP26f7YmCnRiirEDF/+yl88MMplJZrpY5GpFcsL0RERqSelQKR43rg3cC2EARgQ2Imxq1OxM2CEqmjEekNywsRkZERBAFhL7bB6pCesLNS4Mjl2xi2NA4nr96ROhqRXuilvERGRqJFixZQKpXw9PREbGzsE/c/cOAAPD09oVQq0bJlS6xYsUIfMYmIjMpLHi7YMcMPLZ3qIVtdjJErDmFb8lWpYxHpnM7Ly5YtWzB79mwsWLAAKSkpCAgIwMCBA5GZWf1As0uXLmHQoEEICAhASkoK3n//fcyaNQvbtm3TdVQiIqPT2tkWP8zwQz8PZ5SWa/HOdyewaFcqyis4DoZMlyDqeMpGLy8v9OjRA8uXL6/c5uHhgaCgIISHhz+y/9/+9jfs3LkTaWlpldtCQ0Nx4sQJHDp06KnfT6PRQKVSQa1Ww97evm5+CCIiA6fVilj8y3ks+fUiAMC3lQOWju2BhvUsJU5GVDO1+fzW6ZWX0tJSJCcnIzAwsMr2wMBAJCQkVHvMoUOHHtl/wIABOHr0KMrKyh7Zv6SkBBqNpsqLiMjcyGQC5ga2w4rxPWBjKUdCej6GRsThzHW11NGI6pxOy0teXh4qKirg4uJSZbuLiwtycnKqPSYnJ6fa/cvLy5GXl/fI/uHh4VCpVJUvNze3uvsBiIiMzMudGmPHdD+4O9jg2p17eHV5AnaduC51LKI6pZcBu4IgVPm3KIqPbHva/tVtB4D58+dDrVZXvrKysuogMRGR8WrXyA47Z/ijT1snFJdpMTMqBf/66SwqtFzYkUyDTsuLo6Mj5HL5I1dZcnNzH7m68lCjRo2q3V+hUMDBweGR/a2srGBvb1/lRURk7lQ2Flg3sRem9W0JAFhxIB2Tvj4CddGjt9+JjI1Oy4ulpSU8PT0RExNTZXtMTAx8fX2rPcbHx+eR/fft24eePXvCwsJCZ1mJiEyNXCZg/kAPLHm9O5QWMhw8fxPDl8Xh/I0CqaMRPRed3zaaO3cuVq9ejbVr1yItLQ1z5sxBZmYmQkNDAdy/7RMSElK5f2hoKK5cuYK5c+ciLS0Na9euxZo1a/Duu+/qOioRkUka1rUJtr3li6b1rXE5vwgjlsVj7+nqxx0SGQOdl5fRo0dj8eLFWLRoEbp164aDBw8iOjoa7u7uAIDs7Owqc760aNEC0dHR+P3339GtWzd88sknWLJkCV599VVdRyUiMlkdm6iwa6Y/fFo6oLC0AqEbkvFFzHloOQ6GjJDO53nRN87zQkT0eGUVWvwzOg3r4i8DAPp5OON/o7vBTsnb8iQtg5nnhYiIDIuFXIaPhnbEf1/rCkuFDL+k5SJoWTzSb96VOhpRjbG8EBGZoZGervhumg8a2SuRfrMQQUvj8evZG1LHIqoRlhciIjPV1a0+ds70Q0/3BigoKceb3xzFst8uwsRGE5AJYnkhIjJjznZKbJrijXFezSCKwH9+PocZm46hsKRc6mhEj8XyQkRk5iwVMvxjRGeEv9IZFnIB0ady8EpkAq7kF0odjahaLC9ERAQAeL13M2ye6g0nOyucu1GAYUvjcfD8TaljET2C5YWIiCp5ujfErjB/dHWrD/W9MkxcdxhfHUznOBgyKCwvRERURSOVElumeuM1T1doReCf0Wfx9ubjuFdaIXU0IgAsL0REVA2lhRyfjeyCRcM7QiETsPPEdby6PAFXbxdJHY2I5YWIiKonCAJCfJpjw2QvONSzRGq2BsOWxiMhPU/qaGTmWF6IiOiJvFs6YOdMf3Rqao9bhaUIXnMYa+MucRwMSYblhYiInqppfWt8H+qLEd2bokIrYtHuVLzz3QkUl3EcDOkfywsREdWI0kKOL0Z1xQeDPSCXCdh+7BpGrTyE63fuSR2NzAzLCxER1ZggCJgc0BLfvtEbDWwscPKqGsOWxuHwpVtSRyMzwvJCRES15tfaETvD/OHR2B55d0sxdlUi1h+6zHEwpBcsL0RE9EzcGtpg21s+GNKlMcq1Iv7+4xnM23YKJeUcB0O6xfJCRETPzMZSgYjXu2PewPaQCcCWo1kYvTIRNzTFUkcjE8byQkREz0UQBIT2bYV1k3rDXqnA8aw7GBIRh+QrHAdDusHyQkREdaJvWyfsmumPdi52uFlQgjFfJWJTUqbUscgEsbwQEVGdcXeoh+3TfTGwUyOUVYh4f8cpvL/jFErLtVJHIxPC8kJERHWqnpUCkeN64L0B7SAIwKakTLy+KhG5HAdDdYTlhYiI6pwgCJjxQmusndALdkoFkq/cxtClcTiWeVvqaGQCWF6IiEhnXmjvjJ1h/mjjbIsbmhKMWZmIrUeypI5FRo7lhYiIdKqFYz3smOGHAR1dUFqhxV+3ncQHP3AcDD07lhciItI5WysFlo/zxDv920IQgA2JmRi3OhG5BRwHQ7XH8kJERHohkwmY+VIbrA7pCTsrBY5cvo1hEfE4nnVH6mhkZFheiIhIr17ycMEPYX5o5VQPOZpijFpxCFuPchwM1RzLCxER6V0rJ1v8MMMP/Ts8GAfz/Ul8+ONplFVwHAw9HcsLERFJwk5pgZXjPTGnX1sAwLeHrmDcqiTcLCiROBkZOpYXIiKSjEwm4O1+bbBmwv1xMIcv38LQiDiOg6EnYnkhIiLJPTIOZiXHwdDjsbwQEZFBqDIOpvz+OJi//3Ca88HQI1heiIjIYPx5HMz6xCucD4YewfJCREQG5c/jYB7OB5PCdZHoAZYXIiIySA/HwbR2tkWOphijVyZiy5FMqWORAWB5ISIig9XKyRY7pvsi8MF8MH/bdgoLdnBdJHPH8kJERAbNTmmBFeP/f12kjUmZeH1VInI1HAdjrlheiIjI4D1cF2nNhJ6wUyqQfOU2hkTEIfkKx8GYI5YXIiIyGi+2d8HOMH+0cbZFbkEJxnx1CJuSOA7G3LC8EBGRUWnhWA87ZvhhYKdGKKsQ8f6OU5i37SRKyiukjkZ6wvJCRERGx9ZKgchxPfDegHYQBGDzkSyMXpmIHDXHwZgDlhciIjJKgiBgxgutsW5iL9grFTiedQdDIuJw+NItqaORjrG8EBGRUftLO2fsmumP9o3skHe3BGNXJeKbhMsQRVHqaKQjLC9ERGT03B3qYft0Xwzt2gTlWhEf7TyDd747geIyjoMxRSwvRERkEmwsFVgyphsWDPKATAC2H7uGkSsScPV2kdTRqI6xvBARkckQBAFT+rTEhje90LCeJU5f02DY0ngkXMyTOhrVIZYXIiIyOb6tHbEzzA+dmtrjVmEpxq9JwqqDGRwHYyJYXoiIyCS5NrDB96G+eLWHK7Qi8I/oNMzafBxFpeVSR6PnxPJCREQmS2khx39f64JFwztCIROw68R1vBKZgCv5hVJHo+fA8kJERCZNEASE+DRH1FRvONpa4WxOAYZGxOG3s7lSR6NnxPJCRERmoVfzhtgzyx89mtWHprgcb3xzBEv2X4BWy3Ewxkan5eX27dsIDg6GSqWCSqVCcHAw7ty589j9y8rK8Le//Q2dO3dGvXr10KRJE4SEhOD69eu6jElERGbCxV6JzVN9MM6rGUQR+CLmPKauT4amuEzqaFQLOi0vY8eOxfHjx7F3717s3bsXx48fR3Bw8GP3LyoqwrFjx/D3v/8dx44dw/bt23H+/HkMGzZMlzGJiMiMWCpk+MeIzvjs1S6wVMjwS9oNDF8aj/M3CqSORjUkiDp6biwtLQ0dOnRAYmIivLy8AACJiYnw8fHB2bNn0a5duxq9z5EjR9C7d29cuXIFzZo1e+r+Go0GKpUKarUa9vb2z/UzEBGRaTt59Q5C1yfjuroYNpZy/GdkVwzu0ljqWGapNp/fOrvycujQIahUqsriAgDe3t5QqVRISEio8fuo1WoIgoD69etX+/WSkhJoNJoqLyIiopro4lofu2b6w7eVA4pKKzBj0zGER6ehvEIrdTR6Ap2Vl5ycHDg7Oz+y3dnZGTk5OTV6j+LiYsybNw9jx459bAsLDw+vHFOjUqng5ub2XLmJiMi8ONha4ds3emNan5YAgJUHMzBh3WHk3y2ROBk9Tq3Ly8KFCyEIwhNfR48eBXD/8bQ/E0Wx2u1/VlZWhjFjxkCr1SIyMvKx+82fPx9qtbrylZWVVdsfiYiIzJxCLsP8QR5YNrYHbCzliL+Yj6ERcTh59Y7U0agaitoeEBYWhjFjxjxxn+bNm+PkyZO4cePGI1+7efMmXFxcnnh8WVkZRo0ahUuXLuHXX3994r0vKysrWFlZ1Sw8ERHREwzu0hhtXGwRuj4ZGXmFGLniED4d3gmjevGqviHR+YDdpKQk9O7dGwCQlJQEb2/vJw7YfVhcLly4gN9++w1OTk61+r4csEtERM9LU1yGuVtO4Je0+/8R/nrvZlg4rAOsFHKJk5kugxiw6+HhgZdffhlTpkxBYmIiEhMTMWXKFAwZMqRKcWnfvj127NgBACgvL8fIkSNx9OhRbNy4ERUVFcjJyUFOTg5KS0t1FZWIiKgKe6UFvgr2xDv920IQgKjDmRi1MhHX79yTOhpBx/O8bNy4EZ07d0ZgYCACAwPRpUsXrF+/vso+586dg1qtBgBcvXoVO3fuxNWrV9GtWzc0bty48lWbJ5SIiIiel0wmYOZLbbBuYi+orC1wIusOhkbEISE9T+poZk9nt42kwttGRERU17JuFWHa+mSkZmsglwn428vtMCWgZY0eQKGaMYjbRkRERKbCraENtk/3xas9XFGhFfHP6LMI25SCuyXlUkczSywvRERENaC0kOO/r3XBJ0GdYCEXsOdUNoKWxeNi7l2po5kdlhciIqIaEgQBwd7u2DzVBy72VriYexdBy+Kx93S21NHMCssLERFRLXm6N8DumQHwatEQd0vKEbrhGP7101kuK6AnLC9ERETPwMnOChsme2GyfwsAwIoD6QhZexh5XFZA51heiIiInpGFXIYPhnTA0rHdYWMpR0L6/WUFUjJvSx3NpLG8EBERPachXZrgxxl+aOlUD9nqYoxemYgNiVdgYrORGAyWFyIiojrQxsUOP87ww8sdG6G0QosPfjiNd787ieKyCqmjmRyWFyIiojpip7TA8vE9MG9ge8gEYNuxq3glMgGZ+UVSRzMpLC9ERER1SBAEhPZthQ2TveBQzxKp2RoMiYjF/geLPNLzY3khIiLSAd9Wjtg9yx/dm9WHprgcb35zFJ/vO4cKLcfBPC+WFyIiIh1prLLGlqk+mODjDgCI+PUiJq47jFuFpRInM24sL0RERDpkqZDh4+GdsHh0N1hbyBF7IQ9DI+JwIuuO1NGMFssLERGRHgR1b4odM3zR3MEG1+7cw2srDmFjEh+nfhYsL0RERHrSvpE9ds70R2AHF5RWaLFgx/3Hqe+V8nHq2mB5ISIi0iN7pQVWBntWeZx6RGQ8LucVSh3NaLC8EBER6dnDx6k3TvaGo60lzuYUYOjSOOw7kyN1NKPA8kJERCQRn1YO2D0zAJ7uDVBQXI6p65Px771cnfppWF6IiIgk1EilxOap3njD7/7q1Mt/T0fwmsO4WcDVqR+H5YWIiEhiFnIZPhx6f3XqepZyHMrIx+AlsThy+ZbU0QwSywsREZGBGNKlCX4M80cbZ1vkFpRgzFeJWB2bwcep/4TlhYiIyIC0drbFDzP8MKxrE1RoRXy6Jw0zNh1DQXGZ1NEMBssLERGRgalnpcCXY7ph0fCOsJALiD6Vg+FL43E2RyN1NIPA8kJERGSABEFAiE9zbJ3mgyYqJTLyChG0LB7bj12VOprkWF6IiIgMWPdmDbB7VgD6tHVCcZkWc7eewPztp1BcZr6z8rK8EBERGbiG9SyxbmIvzO7XBoIARB3OxMgVCcjML5I6miRYXoiIiIyAXCZgdr+2+HpSbzSwscDpaxoMiYhFTOoNqaPpHcsLERGREenb1gl7ZgWge7P60BSXY8q3R/Gvn8xrVl6WFyIiIiPTpL41tkz1wSS/5gCAFQfSMXZ1EnI1xdIG0xOWFyIiIiNkqZDho6EdsWxsD9haKXD40i0MWhKLhIt5UkfTOZYXIiIiIza4S2PsDPND+0Z2yLtbivFrkhCx/wK0WtOdlZflhYiIyMi1dLLFjul+eM3TFVoR+DzmPCZ9fQS3CkuljqYTLC9EREQmwNpSjv+81hWfvdoFVgoZDpy/iSFLYnEs87bU0eocywsREZEJGdXLDTum+6GFYz1cVxdj1IpDWBN3yaQWd2R5ISIiMjEdmthjZ5gfBndujHKtiE92p+KtDcegMZHFHVleiIiITJCd0gJLx3bHx8PuL+6490wOhkbE4fQ1tdTRnhvLCxERkYkSBAETfJvju1BfNK1vjSv5RXhleQI2JF4x6ttILC9EREQmrptbfeyZ5Y9+Hs4oLdfigx9O4+3Nx3G3pFzqaM+E5YWIiMgM1LexxKqQnlgwyANymYCdJ65j2NI4nM3RSB2t1lheiIiIzIQgCJjSpyW2TPVGI3slMm4WImhZPLYeyTKq20gsL0RERGamZ/OGiH47AH3bOqG4TIu/bjuJd747gaJS47iNxPJCRERkhhrWs8S6ib3w3oB2kAnA9mPXMGxpPM7fKJA62lOxvBAREZkpmUzAjBdaI2qKN5ztrHAx9y6GLY3Dd0ezpI72RCwvREREZs6rpQOi3w5AQBtHFJdp8d73J/HOVsO9jcTyQkRERHC0tcI3k3rj3cC2kAnAtmNXDfY2EssLERERAbh/GynsxTbY9KfbSFuPGtbTSCwvREREVIX3n24j/fXBbaRCA5nUjuWFiIiIHvHwNlLl00gp1wxmUjuWFyIiIqrWw6eRNk/1QSN7JdJvFmL40nhEHc6U9DYSywsRERE9Ue8W9ye1+0s7J5Q8WBspI69QsjwKyb4zERERGY2G9SyxdkIvfBWbAZkAtHKylSyLTq+83L59G8HBwVCpVFCpVAgODsadO3dqfPy0adMgCAIWL16ss4xERERUMzKZgNC+rTC1Tytpc+jyzceOHYvjx49j79692Lt3L44fP47g4OAaHfvDDz8gKSkJTZo00WVEIiIiMjI6u22UlpaGvXv3IjExEV5eXgCAVatWwcfHB+fOnUO7du0ee+y1a9cQFhaGn3/+GYMHD9ZVRCIiIjJCOrvycujQIahUqsriAgDe3t5QqVRISEh47HFarRbBwcF477330LFjx6d+n5KSEmg0miovIiIiMl06Ky85OTlwdnZ+ZLuzszNycnIee9y///1vKBQKzJo1q0bfJzw8vHJMjUqlgpub2zNnJiIiIsNX6/KycOFCCILwxNfRo0cBAIIgPHK8KIrVbgeA5ORkfPnll/j6668fu8+fzZ8/H2q1uvKVlWXYK2ESERHR86n1mJewsDCMGTPmifs0b94cJ0+exI0bNx752s2bN+Hi4lLtcbGxscjNzUWzZs0qt1VUVOCdd97B4sWLcfny5UeOsbKygpWVVe1+CCIiIjJatS4vjo6OcHR0fOp+Pj4+UKvVOHz4MHr37g0ASEpKglqthq+vb7XHBAcHo1+/flW2DRgwAMHBwZg0aVJtoxIREZEJ0tnTRh4eHnj55ZcxZcoUrFy5EgAwdepUDBkypMqTRu3bt0d4eDhGjBgBBwcHODg4VHkfCwsLNGrU6IlPJxEREZH50Ok8Lxs3bkTnzp0RGBiIwMBAdOnSBevXr6+yz7lz56BWq3UZg4iIiEyIIEq5spIOaDQaqFQqqNVq2NvbSx2HiIiIaqA2n99cmJGIiIiMCssLERERGRWWFyIiIjIqOnvaSCoPh/BwmQAiIiLj8fBzuyZDcU2uvBQUFAAAlwkgIiIyQgUFBVCpVE/cx+SeNtJqtbh+/Trs7OxqvMSAPmk0Gri5uSErK4tPQ9UAz1ft8ZzVDs9X7fGc1R7P2dOJooiCggI0adIEMtmTR7WY3JUXmUwGV1dXqWM8lb29PX+Ba4Hnq/Z4zmqH56v2eM5qj+fsyZ52xeUhDtglIiIio8LyQkREREaF5UXPrKys8NFHH3El7Bri+ao9nrPa4fmqPZ6z2uM5q1smN2CXiIiITBuvvBAREZFRYXkhIiIio8LyQkREREaF5YWIiIiMCstLHQoPD0evXr1gZ2cHZ2dnBAUF4dy5czU+Pj4+HgqFAt26ddNdSAPyrOerpKQECxYsgLu7O6ysrNCqVSusXbtWD4ml96znbOPGjejatStsbGzQuHFjTJo0Cfn5+XpILK3ly5ejS5culROD+fj44KeffnriMQcOHICnpyeUSiVatmyJFStW6CmtYajtOdu+fTv69+8PJyenyv1//vlnPSaW3rP8nj1kbn/36wrLSx06cOAAZsyYgcTERMTExKC8vByBgYEoLCx86rFqtRohISF46aWX9JDUMDzr+Ro1ahT279+PNWvW4Ny5c4iKikL79u31lFpaz3LO4uLiEBISgjfffBNnzpzBd999hyNHjmDy5Ml6TC4NV1dX/Otf/8LRo0dx9OhRvPjiixg+fDjOnDlT7f6XLl3CoEGDEBAQgJSUFLz//vuYNWsWtm3bpufk0qntOTt48CD69++P6OhoJCcn44UXXsDQoUORkpKi5+TSqe05e8gc/+7XGZF0Jjc3VwQgHjhw4Kn7jh49Wvzggw/Ejz76SOzatavuwxmgmpyvn376SVSpVGJ+fr4ekxmumpyz//znP2LLli2rbFuyZIno6uqq63gGqUGDBuLq1aur/dpf//pXsX379lW2TZs2TfT29tZHNIP1pHNWnQ4dOogff/yxDhMZvpqcM/7df3a88qJDarUaANCwYcMn7rdu3Tqkp6fjo48+0kcsg1WT87Vz50707NkTn332GZo2bYq2bdvi3Xffxb179/QV06DU5Jz5+vri6tWriI6OhiiKuHHjBr7//nsMHjxYXzENQkVFBTZv3ozCwkL4+PhUu8+hQ4cQGBhYZduAAQNw9OhRlJWV6SOmQanJOfszrVaLgoKCp/7dM1U1PWf8u/98TG5hRkMhiiLmzp0Lf39/dOrU6bH7XbhwAfPmzUNsbCwUCvP9n6Om5ysjIwNxcXFQKpXYsWMH8vLyMH36dNy6dctsxr08VNNz5uvri40bN2L06NEoLi5GeXk5hg0bhoiICD2mlc6pU6fg4+OD4uJi2NraYseOHejQoUO1++bk5MDFxaXKNhcXF5SXlyMvLw+NGzfWR2TJ1eac/dnnn3+OwsJCjBo1SscpDUttzhn/7j8/XnnRkbCwMJw8eRJRUVGP3aeiogJjx47Fxx9/jLZt2+oxneGpyfkC7v9XnSAI2LhxI3r37o1Bgwbhiy++wNdff212V19qes5SU1Mxa9YsfPjhh0hOTsbevXtx6dIlhIaG6imptNq1a4fjx48jMTERb731FiZMmIDU1NTH7i8IQpV/iw8mIf/zdlNW23P2UFRUFBYuXIgtW7bA2dlZD0kNR03PGf/u1xFJb1qZqLCwMNHV1VXMyMh44n63b98WAYhyubzyJQhC5bb9+/frKbG0anq+RFEUQ0JCxFatWlXZlpqaKgIQz58/r6uIBqc252z8+PHiyJEjq2yLjY0VAYjXr1/XVUSD9dJLL4lTp06t9msBAQHirFmzqmzbvn27qFAoxNLSUn3EM0hPOmcPbd68WbS2thZ3796tp1SG7XHnjH/36wavV9UhURQxc+ZM7NixA7///jtatGjxxP3t7e1x6tSpKtsiIyPx66+/4vvvv3/q8cautucLAPz8/PDdd9/h7t27sLW1BQCcP38eMpkMrq6uuo4suWc5Z0VFRY9cmpbL5ZXvZ25EUURJSUm1X/Px8cGuXbuqbNu3bx969uwJCwsLfcQzSE86Z8D9Ky5vvPEGoqKizG4s1eM87pyZ+9/9OiNhcTI5b731lqhSqcTff/9dzM7OrnwVFRVV7jNv3jwxODj4se9hTqPOn+V8FRQUiK6uruLIkSPFM2fOiAcOHBDbtGkjTp48WYofQe+e5ZytW7dOVCgUYmRkpJieni7GxcWJPXv2FHv37i3Fj6BX8+fPFw8ePCheunRJPHnypPj++++LMplM3LdvnyiKj56rjIwM0cbGRpwzZ46YmpoqrlmzRrSwsBC///57qX4EvavtOdu0aZOoUCjEZcuWVfmdvHPnjlQ/gt7V9pz9mTn93a8rLC91CEC1r3Xr1lXuM2HCBLFv376PfQ9z+iV+1vOVlpYm9uvXT7S2thZdXV3FuXPnVvnwNmXPes6WLFkidujQQbS2thYbN24sjhs3Trx69ap+w0vgjTfeEN3d3UVLS0vRyclJfOmllyo/UESx+nP1+++/i927dxctLS3F5s2bi8uXL9dzamnV9pz17du32t/JCRMm6D+8RJ7l9+yPzOnvfl0RRNEMrxsTERGR0eLTRkRERGRUWF6IiIjIqLC8EBERkVFheSEiIiKjwvJCRERERoXlhYiIiIwKywsREREZFZYXIiIiMiosL0RERGRUWF6IiIjIqLC8EBERkVFheSEiIiKj8n+2d58TBRej9AAAAABJRU5ErkJggg==", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "a = start[1]\n", + "b = end[1]\n", + "domain = np.linspace(a,b, 100)\n", + "plt.plot(domain, np.sin(domain))" + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "b544a0bc", + "metadata": { + "scrolled": true + }, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Integral in thrid domain: 0.2167727513247394\n" + ] + } + ], + "source": [ + "print('Integral in thrid domain: {}'.format(sin_integral(a,b)))" + ] + }, + { + "cell_type": "markdown", + "id": "faee288f-e3a6-4c01-9458-d7c7d0751cdf", + "metadata": {}, + "source": [ + "## 2. BTC Workflow\n", + "\n", + "Now we present the complete Workflow for the **BTC** for the **AE kernel**." + ] + }, + { + "cell_type": "markdown", + "id": "4725ec35", + "metadata": {}, + "source": [ + "### 2.1. Domain Discretization\n", + "\n", + "The first thing to do for computing the integral in a computer is the discretization of the domain. In the benchmark we always discretize the domain in $2^n$ intervals, with $n \\in \\mathbf{N}$:\n", + "\n", + "$$\\{[x_0, x_1], [x_1, x_2], ..., [x_{2^n-1}, x_{2^n}]\\}$$ \n", + "\n", + "Where\n", + "\n", + "1. $x_{i+1} < x_{i}$\n", + "2. $a = x_0$\n", + "3. $b = x_{2^n}$" + ] + }, + { + "cell_type": "code", + "execution_count": 8, + "id": "c4d9ce5f", + "metadata": {}, + "outputs": [], + "source": [ + "#First integration domain\n", + "a = start[0]\n", + "b = end[0]\n", + "#For fix the number of discretization intervals\n", + "n=4\n", + "domain_x = np.linspace(a, b, 2**n+1)" + ] + }, + { + "cell_type": "markdown", + "id": "8f43650e", + "metadata": {}, + "source": [ + "### 2.2. Function discretization\n", + "\n", + "Now using the domain discretization we need to discretized $f(x)$. In our procedure following arrays should be constructed:\n", + "\n", + "1. $\\Delta x_i = x_{i+1} - x_{i} = \\frac{b-a}{2^n}$\n", + "2. $f_{x_i} = \\frac{f(x_{i+1}) + f(x_{i})}{2}$\n", + "3. $f_{x_i} \\Delta x_i = f_{x_i} \\frac{b-a}{2^n}$\n", + "\n", + "Using these computed arrays the desired integral can be approximated by Riemann sum:\n", + "\n", + "$$S_{[a,b]} = \\sum_{i=0}^{2^n-1} f_{x_i} \\Delta x_i$$\n", + "\n", + "When $\\Delta x_i \\rightarrow 0$ then $I =\\int_a^{b}\\sin(x)dx \\approx S_{[a,b]}$.\n", + "\n", + "Using $\\Delta x_i = \\frac{b-a}{2^n}$ then we can write down:\n", + "\n", + "\n", + "$$S_{[a,b]} = \\sum_{i=0}^{2^n-1} f_{x_i} \\frac{b-a}{2^n} = \\frac{b-a}{2^n} \\sum_{i=0}^{2^n-1} f_{x_i}$$" + ] + }, + { + "cell_type": "code", + "execution_count": 9, + "id": "a30668da", + "metadata": {}, + "outputs": [], + "source": [ + "#The selected fucntion\n", + "f = np.sin\n", + "\n", + "#Discretization of the selected function\n", + "f_x = []\n", + "x_ = []\n", + "for i in range(1, len(domain_x)):\n", + " step_f = (f(domain_x[i]) + f(domain_x[i-1]))/2.0\n", + " #print(i)\n", + " f_x.append(step_f)\n", + " x_.append((domain_x[i] + domain_x[i-1])/2.0)\n", + "f_x = np.array(f_x)\n", + "x_ = np.array(x_)\n" + ] + }, + { + "cell_type": "code", + "execution_count": 10, + "id": "39cabf6d", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "" + ] + }, + "execution_count": 10, + "metadata": {}, + "output_type": "execute_result" + }, + { + "data": { + "image/png": "iVBORw0KGgoAAAANSUhEUgAAAiMAAAGdCAYAAADAAnMpAAAAOXRFWHRTb2Z0d2FyZQBNYXRwbG90bGliIHZlcnNpb24zLjguMiwgaHR0cHM6Ly9tYXRwbG90bGliLm9yZy8g+/7EAAAACXBIWXMAAA9hAAAPYQGoP6dpAABIAklEQVR4nO3de1xUdf4/8NeZ4TIqzBgoN0XCu0TeIExcvl1VzNXcra+aKVnulrpd1LXMr1um9Vu0Xc2yxC5eMs2062LrkuxuKV6KRGylcbOUQGUAEZsBkcvMOb8/RtBhBpzBmTlzeT0fDx4+OHxmeHPUmRef8znvjyBJkgQiIiIimSjkLoCIiIj8G8MIERERyYphhIiIiGTFMEJERESyYhghIiIiWTGMEBERkawYRoiIiEhWDCNEREQkqwC5C7CHKIooKytDaGgoBEGQuxwiIiKygyRJqKmpQUxMDBSKtuc/vCKMlJWVITY2Vu4yiIiIqANOnz6Nnj17tvl1rwgjoaGhAMw/jFqtlrkaIiIisofBYEBsbGzL+3hbvCKMNF+aUavVDCNERERe5lpLLLiAlYiIiGTFMEJERESyYhghIiIiWXnFmhF7SJIEo9EIk8kkdyl+T6lUIiAggLdhExGRXXwijDQ2NkKn06Gurk7uUuiyzp07Izo6GkFBQXKXQkREHs7rw4goiiguLoZSqURMTAyCgoL4G7mMJElCY2Mjzp07h+LiYvTr16/dRjdEREReH0YaGxshiiJiY2PRuXNnucshAJ06dUJgYCBKSkrQ2NgIlUold0lEROTBfOZXVv727Vn490FERPby+pkRIiIi6hiTKCG/uBqVNfWICFUhJT4MSoX7lzowjHiZF154AZ999hmOHj1q92Nuv/12DB06FGvWrJG1DiIi8hw5RTos26WFTl/fcixao8LSCQlIT4x2ay0MI5d5Sjq8loULF+KJJ55w6DGffPIJAgMDXVQRERF5m5wiHeZsPQKp1fFyfT3mbD2CrOnD3RpIGEbgWemwLZIkwWQyISQkBCEhIQ49NiwszEVVERGRtzGJEpbt0loFEQCQAAgAlu3SYnRClNt+Kff7VYbN6fDqIAJcSYc5RTqXfe+GhgY8+eSTiIiIgEqlwq9+9St8++23AICvvvoKgiDgiy++QHJyMoKDg5GXl4cXXngBQ4cObXkOo9GIJ598El27dkV4eDgWLVqEhx56CJMmTWoZc/vtt2PevHktn994443485//jEceeQShoaHo1asX3nrrLYvaFi1ahP79+6Nz587o3bs3nnvuOTQ1NbnsXBARkXvkF1dbveddTQKg09cjv7jabTX5XBiRJAl1jUa7Pmrqm7A0+/s20yEAvJCtRU19k13PJ0m2nqltzzzzDD7++GO8++67OHLkCPr27YuxY8eiurraYkxmZiaOHz+OwYMHWz3HypUrsW3bNmzatAkHDhyAwWDAZ599ds3vvWrVKiQnJ6OwsBBz587FnDlz8N///rfl66Ghodi8eTO0Wi1effVVvP3223jllVcc+vmIiMjzVNa0HUQ6Ms4ZfO4yzaUmExKe/8IpzyUBKDfU4+YX9tg1Xrt8LDoH2XdKL168iKysLGzevBnjxo0DALz99tvIzc3Fhg0bcMsttwAAli9fjtGjR7f5PGvXrsXixYvxm9/8BgDw+uuvY/fu3df8/vfccw/mzp0LwDwL8sorr+Crr77CwIEDAQB/+tOfWsbeeOON+OMf/4gdO3bgmWeesevnIyIiNxJNQMlBoLYCCIkE4lIBhdJiSH2TCV98X471X5206ykjQt3XI8rnwoi3OHnyJJqamjBq1KiWY4GBgUhJScHx48dbwkhycnKbz6HX61FRUYGUlJSWY0qlEklJSRBFsd3vf/UsiyAIiIqKQmVlZcuxjz76CGvWrMFPP/2E2tpaGI1GqNVqh39OIiJyMW02kLMIMJRdOaaOAdJXAgkTcaKiBtvzS/Fp4Vn8Unfty+0CgCiN+UYOd/G5MNIpUAnt8rF2jc0vrsbMTd9ec9zmh2+x6y+lU6DymmOaNV/Sad26XpIki2NdunS55nPZeo5raX13jSAILQHm66+/xtSpU7Fs2TKMHTsWGo0GH3zwAVatWnXN5yUiIjfSZgM7M4BWCw4kgw7YmYG/aJZgXUVCy/EYjQqTb4lFRGgwlnxaZB571eOa302WTkhw6x2lPhdGBEGw+1JJWr/uiNaoUK6vt7lupDkdpvXr7vS/lL59+yIoKAj79+/HtGnTAABNTU04fPiwxWLT9mg0GkRGRiI/Px9paWkAAJPJhMLCQotFro46cOAA4uLisGTJkpZjJSUlHX4+IiJyAdFknhGx8Q4mQIIoAdN/ycLbitdwx8AoPDCiF/7nqvezsC5BVneSRrHPiPspFQKWTkjAnK1HIMC96bBLly6YM2cOnn76aYSFhaFXr154+eWXUVdXh1mzZuG7776z63meeOIJZGZmom/fvhg4cCDWrl2LCxcuXNdmgX379kVpaSk++OAD3HLLLfj73/+OTz/9tMPPR0RELlBy0PLSTCsKAYjBeeQ/2Ak33GR9yT89MRqjE6I8oseWX4cRwPyXkTV9uCzpcMWKFRBFETNmzEBNTQ2Sk5PxxRdf4IYbbrD7ORYtWoTy8nJkZGRAqVTi0UcfxdixY6FU2n/JqLV7770X8+fPx+OPP46GhgaMHz8ezz33HF544YUOPycRETmPJEkoLTmFODvG3iBeaPNrSoWAkX3CnVdYBwmSo/ejysBgMECj0UCv11stoqyvr0dxcTHi4+Ova3dYb+nAei2iKGLQoEGYPHkyXnzxRdnqcNbfCxGRP7nWe5H+UhP+dvQstuefhqbia3wQ9NK1n/Shz4H4NBdW3bb23r+v5vczI808JR06qqSkBHv27MFtt92GhoYGvP766yguLm5Zh0JERN6hrW7gz/86Ad1Dg7E9/zT+fqwM9U3mmw1UAQm4ENAdXY1VENpa+aiOMd/m6+EYRrycQqHA5s2bsXDhQkiShMTERPzzn//EoEGD5C6NiIjs1NZeMTp9PeZsO2JxrH9kCB5I6YXfDOuBrj+vvnw3TRsrH9NXWPUb8UQMI14uNjYWBw4ckLsMIiLqoPb2irna/cN74IERcRjeq+uVmxQSJgKTt7TRZ2SF+etegGGEiIhIRtfaK6bZfUmxSIqzcYNDwkRg4PhrdmD1ZAwjREREMqlvMuGTwjN2jW13rxiFUrZFqs7AMEJERORs19gr5sLFRmw5VIIth37G+YuNdj2lO/eKcTeGESIiImdqZ6+Y01F34528U9hx+HTLXTExGhVqG4yoqTe22w3cnXvFuBvDCBERkbO0u1fMDPy5aR7+YTJvbnpTjBqP3dYH9yRG4Z/HK2TpBu4pGEaIiIicwY69Yp4LeA91vdPx6G39kNonvOWuGDm7gXsChdwF+Kvbb7/d7g3xbNm8eTO6du1qceytt95CbGwsFAoF1qxZc131ERGRg+zZK0Y4j3fvNGJU325We4ilJ0Zj/6I7sf33t+LVqUOx/fe3Yv+iO30+iACcGbniGouNPJ3BYMDjjz+O1atX47777oNGo5G7JCIiv3Lpwll0smdgbUWbX/LWbuDXi2EEaHexkbc0jCktLUVTUxPGjx+P6GjfT9FERJ6iXF+PTQeKceKbMmyyZ1lHSKTLa/I2vEzTvNio9dSaQWc+rs122bcWRRHPPPMMwsLCEBUVZbEr7urVq3HzzTejS5cuiI2Nxdy5c1FbW2vzeTZv3oybb74ZANC7d28IgoCff/7ZZXUTEfkLkyjh0Mnz+NvRszh08jxM4pX1ICcqarDww++Q9vK/8ea+U9jb0A+VQjgktJVIBEDdwyv2inE3/54ZaWexkfmYAOQ8a+5s54JLNu+++y4WLFiAb775BocOHcLMmTMxatQojB49GgqFAq+99hpuvPFGFBcXY+7cuXjmmWewbt06q+eZMmUKYmNjcffddyM/Px+xsbHo3r270+slIvIntjaui9KoMC2lFwpLL+DLH861HE+5MQyP/k9vdJNegfDhQ/D2vWLczb/DyDUWGwESYDhrHueCznaDBw/G0qVLAQD9+vXD66+/jn/9618YPXq0xeLW+Ph4vPjii5gzZ47NMNKpUyeEh5uvMXbv3h1RUVFOr5WIyJ+0tXFdub4eq3NPAAAEARibEIVHb+uN4b2a27TfCwjev1eMu/l3GGlnEVGHxjlo8ODBFp9HR0ejsrISAPDll1/iz3/+M7RaLQwGA4xGI+rr63Hx4kV06dLFJfUQEZF9G9d1DlIi+/FfoW9EiPUXfWCvGHfz7zUj9i4ictFio8DAQIvPBUGAKIooKSnBPffcg8TERHz88ccoKCjAG2+8AQBoampySS1ERGRmz8Z1dY0mnKtpaHtA814xN99v/pNBpF3+PTMSl2qeOjPoYHvdiGD+upsXGx0+fBhGoxGrVq2CQmHOizt37nRrDURE/kgUJeQU6ewa2+7GdeQQ/54ZUSjNt+8CgNXqZ/kWG/Xp0wdGoxFr167FqVOn8N5772H9+vVurYGIyJ9IkoQ935fjntfy8O6hErse48sb17mbf4cRwHxtb/IWQN2qN4c6xnxchsVGQ4cOxerVq7Fy5UokJiZi27ZtyMzMdHsdREQ+QzQBxXnAsY/Mf4omAOYQkvfjOUxadxCPvleA/5bXICRIiZDggPZu0EW0j29c526CJEntrdHxCAaDARqNBnq9Hmq12uJr9fX1KC4uRnx8PFSq60ipXt6B1dM47e+FiOh6tdHY8qek5/B//41HfnE1AKBToBIPj7oRj/5Pb3x96jzmbD0CwPbGdVnTh/tFm/br1d7799X8e83I1ZoXGxERke9oYxdd0VCG3v+egxua5iEo4FZMHxGHObf3QffQYADcuM7dGEaIiMg3tdPYUgFABLCyy/uom/M0YsKsb9FNT4zG6IQo5BdXo7KmHhGh5kszSoU9Pd/JEQwjRETkm+zYRbdrUyW66guBMNsz4/66cZ27cQErERH5pOqK0/YNdFFjS7IfZ0aIiMinVNbUY92XJ/FTfhm22vMux110ZeczYcQLbgryK/z7ICJXMIlSm2s4LlxsxPp9J/HuwZ9R3yRCgQE4H9QNYeJ5CB7U2JKseX0YaW6pXldXh06dOslcDTWrq6sDYN3ynoioo2ztohutUeGZsQNQUl2HDXnFqGkwAgCG9eqKp8cMQHjjK5fvpuEuup7M68OIUqlE165dWzaY69y5MwSBK53lIkkS6urqUFlZia5du0Kp5H9yIrp+be2iq9PXY/7O71o+T4hWY+HY/rhjQMTl94LLjS25i65H8/owAgBRUVEA0BJISH5du3Zt+XshIroe9uyiq1QIWDN5KMYPjoai9a233EXX4/lEGBEEAdHR0YiIiOCuth4gMDCQMyJE5DT27KJrEiV0Cw22DiLN2NjSo/lEGGmmVCr5JkhE5GPs3R2Xu+h6L/YZISIij1VpqMfHR87YNZa76HqvDoWRdevWtWyAlpSUhLy8vHbHb9u2DUOGDEHnzp0RHR2Nhx9+GOfPn+9QwURE5Pvqm0x448ufcMdfv8K+E1XtjuUuut7P4TCyY8cOzJs3D0uWLEFhYSHS0tIwbtw4lJaW2hy/f/9+ZGRkYNasWfj+++/x4Ycf4ttvv8Xvfve76y6eiIi8lGgCivOAYx+Z/xRNAMx35OUU6TD6lb34yxc/4GKjCUNju+KZsQMg4Mquuc2aP186IYF7xngxQXKwO9WIESMwfPhwZGVltRwbNGgQJk2ahMzMTKvxf/3rX5GVlYWTJ0+2HFu7di1efvllnD5tX6tee7cgJiIiL6DNtnmr7ekRL+Dp73vh61PVAIBIdTCeHTcQ9w7pAYVCaLPPCHfR9Vz2vn87tIC1sbERBQUFePbZZy2OjxkzBgcPHrT5mNTUVCxZsgS7d+/GuHHjUFlZiY8++gjjx4935FsTEZEv0GZfbkJm+XuwZChDjz2PQtM0D8EBt+Kx/+mN2bf3QeegK29T3EXXdzkURqqqqmAymRAZadnHPzIyEuXl5TYfk5qaim3btmHKlCmor6+H0WjExIkTsXbt2ja/T0NDAxoaGlo+NxgMjpRJRESeSDSZZ0RsdAxp7o+a2Xkb6uYsRM/wUJtPwV10fVOHFrC27nAqSVKbXU+1Wi2efPJJPP/88ygoKEBOTg6Ki4sxe/bsNp8/MzMTGo2m5SM2NrYjZRIRkScpOWh5aaYVhQCEGc+hp+Go+2oij+BQGOnWrRuUSqXVLEhlZaXVbEmzzMxMjBo1Ck8//TQGDx6MsWPHYt26ddi4cSN0Op3NxyxevBh6vb7lw961JURE5MFqK5w7jnyGQ2EkKCgISUlJyM3NtTiem5uL1FTbux7W1dVBobD8Ns2NydpaOxscHAy1Wm3xQURE3uuXukZs/s8l+waH2P7llnyXwx1YFyxYgBkzZiA5ORkjR47EW2+9hdLS0pbLLosXL8bZs2exZcsWAMCECRPw+9//HllZWRg7dix0Oh3mzZuHlJQUxMTEOPenISIij2I0iXg/vxSrc0/AUNcVY4LDECVUt/GbsGDewC7O9i+35LscDiNTpkzB+fPnsXz5cuh0OiQmJmL37t2Ii4sDAOh0OoueIzNnzkRNTQ1ef/11/PGPf0TXrl1x5513YuXKlc77KYiISBYmUWrz7pa8H8/hxc+1OFFRCwAYEKmBfshLiNn3h8uPvnp2/PK6w/QV3MDODzncZ0QO7DNCROR52ur7Mfu2Psj7sQr/PG5e+3FD50AsGN0fD6T0QoBS0UafkR7mIJIw0d0/BrmQve/fDCNEROSwnCId5mw9YuMm3SuUCgEZI+Mw767+0HQOtPyiaDLfXVNbYV4jEpfKGREf5JKmZ0RERCZRwrJd2naDSHCAAtmPj8KAqDbegBRKID7NJfWR9+GuvURE5JD84mqLSzO2NBhFVF9sclNF5O0YRoiIyCGVNe0HEUfHETGMEBGR3SRJws9VF+0aGxGqcnE15Cu4ZoSIiOxS9sslPP+3IvzzeGW74wQAURrzbb5E9mAYISKidplECZsP/oxVe35AXaMJgUoBdw+KRE6ReWsQG91CsHRCAnfTJbsxjBAR+bt2brMtOqvH4k+O4dhZPQAgOe4GZP72ZvSLDLXZZyRKo8LSCQlIT4yW5Uch78QwQkTkz2w2IItB/d1/xl9LB2DjgWKIEhCqCsDicYMw9ZZYKC7PeKQnRmN0QlSbHViJ7MUwQkTkr7TZwM4MoFXHEMmgQ/AnM3G6cR5EKQXjB0dj6a8TEKG2XpCqVAgY2SfcTQWTr2IYISLyR6LJPCNio3WZAAmiBCwP2oqpk2fjjgReciHXYhghIvJHJQctL820ohCASFQhstNPABhGyLXYZ4SIyB/VVjh3HNF1YBghIvJDjZ262zcwJNK1hRCBYYSIyO8c/KkK6Z80oUwKg9jmbncCoO5hvs2XyMUYRoiI/ET1xUb8ced3mPbONzhV3YBXA2dBEARIaH0r7uXP01e09BshciUuYCUi8hEmUbLZ80OSJHxy5Cxe+rsWF+qaIAjAjFvjsHDsGAinBtvsM4L0FUDCRPl+GPIrDCNERD7AVjfUaI0Ks2/rgz3achz46TwAYEBkKDLvuxnDe91gHpQwERg4vs0OrETuIEiS1OYVQ09hMBig0Wig1+uhVqvlLoeIyKPkFOkwZ+sRGx1DrggOUOCpu/vh92m9EajkFXpyD3vfvzkzQkTkxUyihGW7tNcMIv94Kg29u4e4rS4iRzAeExF5sfziaotLM7Y0GEVUGBrcVBGR4xhGiIi8WGVN+0HE0XFEcmAYISLyYl07B9o1LiLUepM7Ik/BNSNERF6q6KweL32ubXeMACBKY77Nl8hTMYwQEXmZJpOIN778Ca//+ycYRQmhqgDU1BshwHIP3uZWZksnJECpaN3YjMhzMIwQEXmRExU1WLDzKIrOGgAA99wchZcm3Yz84vNWfUaiNCosnZCA9ETuukuejWGEiMjTiCarJmQmKPBO3ims2nMCjSYRmk6BWH7vTZg4JAaCICA9MRqjE6JsdmAl8nQMI0REnkSbbdWe3RgSjVeUs/BGRQIA4M6BEcj87c2IVFsuSlUqBIzsE+7WcomcgWGEiMhTaLOBnRlAqxZmihod/oiXcDb4j0j99cP43+SeEATOeJDv4K29RESeQDSZZ0Rs9FJVCAAE4K+h2zE5KYZBhHwOwwgRkScoOWi5c24rCgABtWXmcUQ+hmGEiMgT1FY4dxyRF2EYISLyAIcq7VzCFxLp2kKIZMAwQkQkowsXG/H4+0fwYK4SZVIYxDZHCoC6BxCX6sbqiNyDYYSISCb/Ol6BMWv24fP/6CAolDg8aBEEXF6tauHy5+krAIXS3WUSuRxv7SUichGTKNlsQmaob8JLn2ux8/AZAEDfiBCsnjwEg3veA2hjrPqMQB1jDiIJE2X6SYhci2GEiMgFcop0Vu3ZozUqTEmOxYcFZ3D2l0sQBOD3ab2xYHR/qAIvz3gkTAQGjrfqwMoZEfJlgiRJ1je1exiDwQCNRgO9Xg+1Wi13OURE7cop0mHO1iM2OoZc0SusM1ZNHoJbbuRuuuS77H3/5swIEZETmUQJy3Zp2w0inYOU+PyJX0HdKdBtdRF5Mi5gJSJyovziaotLM7bUNZrwfZnBTRUReT6GESIiJ6qsaT+IODqOyB8wjBAROVFEqOragxwYR+QPuGaEiMhJLjYY8VHB6XbHCACiNObbfInIjGGEiMgJis7q8cT2QhRXXYQA8967zX82a25ltnRCApQK7rxL1IxhhIjoOoiihI0HirEy579oMkmI0aiwZuowVF9ssOozEqVRYemEBKQnRstYMZHnYRghIroW0WSzCVlVbQMWfvgdvvrhHABg7E2RWHnfYHTtHAQAGJ0QZbMDKxFZYhghImqPNttme3btkCXIOBSFqtoGBAco8NyvE/DgiF4QhCthQ6kQMLJPuAxFE3kXhhEiorZos4GdGUCrFmaSoQwD9/0BSU3z8HPkXXjtgWEYEBUqT41EPoBhhIjIFtFknhGx0Uu1eWHqyi7vQzV3CVTBQe6ujsinsM8IEZEtJQctL820ohCArk2VUJV948aiiHwTwwgRkS21Fc4dR0RtYhghIrIlJNK544ioTQwjREStiKKEDaejoJPCILa5/a4AqHuYb/MlouvCMEJEdJVzNQ14ePO3eHH3CbzQlAFBACS07g1y+fP0FYBC6fYaiXwNwwgR0WX7TpzDuFfzsPfEOQQHKPA/9z4CTN4CQd2qY6o6Bpi8BUiYKE+hRD6Gt/YSkd9rNIpYtecHvLnvFABgQGQo1k4bhv6RoQDigIG/ttmBlYicg2GEiPyGSZSs2rOfrq7Dkx8U4j9n9ACAGbfGYcn4QVAFXhU2FEogPk2mqol8H8MIEfmFnCKd1cZ1XTsF4lKTCQ1GEZpOgXj5/sEYe1OUjFUS+SeGESLyeTlFOszZesSql+ovl5oAAH0jQrDlkRTEdO3k/uKIiAtYici3mUQJy3ZpbTR1v+JigxGRapXbaiIiSwwjROTT8ourLS7N2KLT1yO/uNpNFRFRawwjROTTKmvaDyKOjiMi5+tQGFm3bh3i4+OhUqmQlJSEvLy8dsc3NDRgyZIliIuLQ3BwMPr06YONGzd2qGAiIkc0mdq7QHNFRCgv0xDJxeEFrDt27MC8efOwbt06jBo1Cm+++SbGjRsHrVaLXr162XzM5MmTUVFRgQ0bNqBv376orKyE0Wi87uKJiNrz+X/K8Pxnx9odIwCI0phv8yUieQiSJNn3a8NlI0aMwPDhw5GVldVybNCgQZg0aRIyMzOtxufk5GDq1Kk4deoUwsI69p/dYDBAo9FAr9dDrVZ36DmIyH80GkVk/uM4Nh34GQDQPzIEJypqIQAWC1mbm7xnTR+O9MRWXVaJ6LrZ+/7t0GWaxsZGFBQUYMyYMRbHx4wZg4MHD9p8THZ2NpKTk/Hyyy+jR48e6N+/PxYuXIhLly61+X0aGhpgMBgsPoiI7KHTX8LUtw61BJE5t/fB7ifTsH76cERpLC/FRGlUDCJEHsChyzRVVVUwmUyIjLTcMjsyMhLl5eU2H3Pq1Cns378fKpUKn376KaqqqjB37lxUV1e3uW4kMzMTy5Ytc6Q0IiLs/7EKT35QiOqLjQhVBWD15KEYnWB+vUpPjMbohCirDqxKRetN8IjI3TrU9EwQLP/zSpJkdayZKIoQBAHbtm2DRqMBAKxevRr3338/3njjDXTqZN1kaPHixViwYEHL5waDAbGxsR0plYh8kWiy2CtGjB2JdfuKsSr3BCQJuClGjawHk9ArvLPFw5QKASP7hMtUNBG1xaEw0q1bNyiVSqtZkMrKSqvZkmbR0dHo0aNHSxABzGtMJEnCmTNn0K9fP6vHBAcHIzg42JHSiMhfaLOBnEWAoazl0C/Kbjh2aTokKQVTb4nFCxNvstxbhog8mkNrRoKCgpCUlITc3FyL47m5uUhNTbX5mFGjRqGsrAy1tbUtx06cOAGFQoGePXt2oGQi8lvabGBnhkUQAYCuxipkBa7BtlEVWHHfYAYRIi/jcJ+RBQsW4J133sHGjRtx/PhxzJ8/H6WlpZg9ezYA8yWWjIyMlvHTpk1DeHg4Hn74YWi1Wuzbtw9PP/00HnnkEZuXaIiIbBJN5hkRG43dFYL58vGoH/9qHkdEXsXhNSNTpkzB+fPnsXz5cuh0OiQmJmL37t2Ii4sDAOh0OpSWlraMDwkJQW5uLp544gkkJycjPDwckydPxksvveS8n4KIfF/JQasZkasJkADDWfO4+DQ3FkZE18vhPiNyYJ8RIsKxj4CPZ1173H0bgJvvd309RHRNLukzQkQkl2/O2TmRG2J7MT0ReS6GESLyaE0mES9+rsUDe5Qok8IgtjlSANQ9gDjbi+mJyHMxjBCRxyrX1+OBt77Ghv3FEKHAof5PQ4CAK43cm13+PH0FoOCdNETepkNNz4iIXO3gySo8ub0QVbXmbqp//d8hGHvTeEDb06rPCNQx5iCSMFG+gomowxhGiEg2JlGyas8uAMjaexKr9vwAUQIGRauxfvpwxIV3MT8oYSIwcLxFB1bEpXJGhMiLMYwQkSxyinRYtksLnb6+5VikOhgRocE4dta8Oeb/JvXEi5MSrZuYKZS8fZfIhzCMEJHb5RTpMGfrEav2ZRWGBlQYGhCgEPD/fpOIKbf0kqU+InIvLmAlIrcyiRKW7dLa6KN6RdfOgbg/iZtjEvkLhhEicqv84mqLSzO2VNU2Ir+42k0VEZHcGEaIyK0qa9oPIo6OIyLvxzBCRG4VEapy6jgi8n5cwEpEbvXTuZp2vy4AiNKYb/MlIv/AMEJEbtFoFLFs1/fY9s2VXb0FwGIha3Nf1aUTEqBUtO6ySkS+imGEiFzuXE0D5m4rwLc/X4AgAE+PHYD48C5Y/rlln5EojQpLJyQgPTFaxmqJyN0YRojIpf5z5hc89l4BdPp6hKoC8NrUYbhjYAQAYMxNUVYdWDkjQuR/GEaIyGU+LTyDZz8+hgajiD7du+CtjGT06R7S8nWlQsDIPuEyVkhEnoBhhIiun2iy2CvG2PNWrPjiR7yzvxgAcNfACLwydSjUqkCZCyUiT8QwQkTXR5tttYuuXtkNpy9NB5CCx+/oiwWj+0PByy9E1AaGESLqOG02sDMDaNXc/QZjFbIC16Bw5GtIGjtentqIyGuw6RkRdYxoMs+I2NhlRiEAgiAgSbvSPI6IqB0MI0TUMSUHLS7NtCZAAgxnzeOIiNrBMEJEHVNb4dxxROS3GEaIqEN0Jo19A0MiXVsIEXk9hhEicthXP1Qi/TMjyqQwiG2OEgB1DyAu1Y2VEZE3YhghIrtJkoT1e0/i4c3fQl8vYmvXuRAg4MquMs0uf56+AlAo3V0mEXkZ3tpLRHa51GjCoo//g+zvzItWH0iJxVMT0yGcGGDVZwTqGHMQSZgoU7VE5E0YRojoms5cqMNj7xXg+zIDAhQClk68CdNH9IIgCObAMXC8RQdWxKVyRoSI7MYwQkQtTKJktXHdtz9XY+62I6i+2IjwLkFY9+BwjOjdaj8ZhRKIT5OnaCLyegwjRAQAyCnSYdkuLXT6+pZjalUAahuMECUgsYcab85IRo+unWSskoh8EcMIESGnSIc5W49Y9VI11BsBAMlxN+C9WSPQKYiXXojI+Xg3DZGfM4kSlu3S2mjqfsXZXy4hKIAvF0TkGnx1IfJz+cXVFpdmbNHp65FfXO2miojI3zCMEPm5ypr2g4ij44iIHMUwQuTnIkKD7RyncnElROSvuICVyI81GE3Y8e3pdscIAKI05tt8iYhcgWGEyE+dr23AY+8V4HDJBSgEQJTMwePqhazNTd6XTkiAUtG65TsRkXMwjBD5oR8ravDIu9/idPUlhKoCsO7B4bjYYLTqMxKlUWHphASkJ0bLWC0R+TqGESI/s/fEOTy+7QhqGozoFdYZG2fegr4RIQCA0QlRVh1YOSNCRK7GMELkR949+DOW7foeogSk3BiG9TOSENYlqOXrSoWAkX3C23kGIiLnYxgh8kWiyWLjOmPPW7F89w/YcqgEAHB/Uk/8v98kIjiAHVWJSH4MI0S+RpsN5CwCDGUth/TKbqi4NB1AChalD8Ts23qbd9wlIvIADCNEvkSbDezMAFo1d7/BWIWswDX4LnUtht3eR57aiIjawKZnRL5CNJlnRGzsMqMQAEEQMOz7FeZxREQehGGEyFeUHLS4NNOaAAkwnDWPIyLyIAwjRL6itsK544iI3IRhhMhHNKi62zcwJNK1hRAROYhhhMgHVBjqMSUHKJPCIFovGblMANQ9gLhUd5ZGRHRNDCNEXq7orB73vn4AR8/WYpXikcu37La+bffy5+krAAV7ixCRZ2EYIfJiX3xfjv9dfwjlhnr06d4FTz6+AMLkLYC61V4y6hhg8hYgYaI8hRIRtYN9Roi8kCRJeGvfKazI+S8kCUjr1w2vTxsOTadAIHwiMHC8RQdWxKVyRoSIPBbDCJGXaTSK+NNnx7Dz8BkAwIxb47B0QgIClFdNdCqUQHyaTBUSETmGYYTIQ5lEyWoHXcOlJszeWoBviquhEIDnf52AmaPi5S6ViOi6MIwQeaCcIh2W7dJCp69vOdY9xLy77rnaRoQEB2DttGG4Y0CEXCUSETkNwwiRh8kp0mHO1iNWTd3P1TYCAMK7BOH939+KAVGh7i+OiMgFeDcNkQcxiRKW7dLa2F3migClgL4RIW6riYjI1RhGiDxIfnG1xaUZWyoMDcgvrnZTRURErscwQuRBKmvaDyKOjiMi8gYMI0QeJCJU5dRxRETegGGEyINEqoOhVLRu5X6FACBaY77Nl4jIVzCMEHmII6UXcP/6QzC1sdNdc0RZOiGh3cBCRORtGEaIPEBOUTkeeOtrVF9sxM09NFh5382I1lheionSqJA1fTjSE6PbeBYiIu/EPiNEMtuwvxgv/V0LSQLuGhiB1x4Yhi7BAbg/KdaqAytnRIjIF3VoZmTdunWIj4+HSqVCUlIS8vLy7HrcgQMHEBAQgKFDh3bk2xL5FHNPke/x4ufmIDL91l54c0YSugSbf0dQKgSM7BOOe4f2wMg+4QwiROSzHA4jO3bswLx587BkyRIUFhYiLS0N48aNQ2lpabuP0+v1yMjIwF133dXhYom8kmgCivOAYx+Z/xRNuNRowtxtBdh04GcAwLPjBuLFexMtN7sjIvITgiRJ7TV7tDJixAgMHz4cWVlZLccGDRqESZMmITMzs83HTZ06Ff369YNSqcRnn32Go0eP2v09DQYDNBoN9Ho91Gq1I+USyUubDeQsAgxlLYdMITH4i+JhrK+8CUFKBVZNHoIJQ2JkLJKIyDXsff926NewxsZGFBQUYMyYMRbHx4wZg4MHD7b5uE2bNuHkyZNYunSpXd+noaEBBoPB4oPI62izgZ0ZFkEEAITaMjyj/3/4reoItv5uBIMIEfk9h8JIVVUVTCYTIiMjLY5HRkaivLzc5mN+/PFHPPvss9i2bRsCAuxbL5uZmQmNRtPyERsb60iZRPITTeYZERu7zCgAQABWhryPlDiNuysjIvI4HbpALQiWC+kkSbI6BgAmkwnTpk3DsmXL0L9/f7uff/HixdDr9S0fp0+f7kiZRPIpOWg1I3I1BYDA2jLzOCIiP+fQrb3dunWDUqm0mgWprKy0mi0BgJqaGhw+fBiFhYV4/PHHAQCiKEKSJAQEBGDPnj248847rR4XHByM4OBgR0oj8iy1Fc4dR0TkwxyaGQkKCkJSUhJyc3Mtjufm5iI1NdVqvFqtxrFjx3D06NGWj9mzZ2PAgAE4evQoRowYcX3VE3mqEOtwfl3jiIh8mMNNzxYsWIAZM2YgOTkZI0eOxFtvvYXS0lLMnj0bgPkSy9mzZ7FlyxYoFAokJiZaPD4iIgIqlcrqOJEvqYtOQb2yG7oaq2C7PYgAqGOAOOsQT0TkbxwOI1OmTMH58+exfPly6HQ6JCYmYvfu3YiLiwMA6HS6a/YcIfJl52oaMOvdbxF9aTqyAtdAggDBYiHr5XSSvgJQKGWpkYjIkzjcZ0QO7DNC3uKnylrM3JSPMxcu4YbOgfjwtir0LXjRcjGruoc5iCRMlK9QIiI3sPf9m3vTEDnJN6fO49H3CqC/1IS48M7Y/HAK4rt1AdKmmO+aqa0wrxGJS+WMCBHRVRhGiJwg+7syLNz5HRpNIob16op3MpIRHnL5jjCFEohPk7dAIiIPxjBCdB0kScL6vaewMue/AICxN0Xi1anDoArkzAcRkb0YRojsZBIl5BdXo7KmHhGhKgzv1RXLPtfi/W/MC7YfGRWPJeMHcXddIiIHMYwQ2SGnSIdlu7TQ6etbjgUHKNBgFCEIwHPjE/DIr+JlrJCIyHsxjBBdQ06RDnO2HrHaZabBKAIAHk3rzSBCRHQdOrQ3DZG/MIkSlu3S2tju7ors78pgEj3+DnkiIo/FMELUjvziaotLM7bo9PXIL652U0VERL6HYYSoHZU17QcRR8cREZE1hhGidkSEqpw6joiIrHEBK1EbJEnCgZNV7Y4RAERpVEiJD3NPUUREPohhhMgGo0nEnz4rwgffnm45JgC2trvD0gkJ7C1CRHQdeJmGqJVLjSbM3lqAD749DYUAvDQpEeunD0eUxvJSTJRGhazpw5GeGC1TpUREvoEzI0RXuXCxEbPe/RZHSn9BUIACr00dhvTEKADA6IQoiw6sKfFhnBEhInIChhGiy85cqMNDG/Nx8txFqFUB2DDzFtxy45W1IEqFgJF9wmWskIjINzGMkP8RTUDJQaC2AgiJBOJS8d/Ki3hoYz4qDA2I1qjw7iMp6B8ZKnelRER+gWGE/Is2G8hZBBjKWg41dI7C+roHUVGfhH4RIXj3kRTEdO0kY5FERP6FC1jJf2izgZ0ZFkEEAAIvlmO1tAp/iNLio9mpDCJERG7GMEL+QTSZZ0Rs7DKjEABBABaKm6BR8b8EEZG78ZWX/EPJQasZkasJAATDWfM4IiJyK4YR8g+1Fc4dR0RETsMwQv4hJNK544iIyGkYRsgvVHdLRpWiG0TrJSOXCYC6BxCX6s6yiIgIDCPkB05X1+H+N7/BkvrpgABIaN019fLn6SsAhdLt9RER+TuGEfJp35fp8dusgzhVdRFF6ttQMfZtCOpWe8moY4DJW4CEifIUSUTk59j0jHzWwZ+q8Oh7BahtMGJgVCjefSQFkWoVMOI+qw6snBEhIpIPwwj5pF3flWHBzqNoMkkYER+GtzKSoekUaP6iQgnEp8lbIBERtWAYIZ+zcX8xln+uBQDcc3MUVk8eClUgZz6IiDwVwwh5LZMoIb+4GpU19YgIVSE57gb8NfcHvLn3FADgoZFxeH7CTVAqWi9YJSIiT8IwQl4pp0iHZbu00OnrW451ClTgUpMIAHgmfQDm3NYHgsAgQkTk6RhGyOvkFOkwZ+sRq11mmoNIxsg4zL29r/sLIyKiDuGtveRVTKKEZbu0Nra7uyJXWwFT293NiIjIwzCMkFfJL662uDRji05fj/ziajdVRERE14thhLxKZU37QcTRcUREJD+GEfIqEaEqp44jIiL5cQEreZW6BmO7XxcARGlUSIkPc09BRER03TgzQl7js8KzeGxrQcvnbWx3h6UTEthbhIjIizCMkFfYdKAY83YchVGU8JthPfD6tGGI0lheionSqJA1fTjSE6PbeBYiIvJEvExDHk2SJKzOPYG1//4JAPDwqBvx3PgEKBQCxiVGW3RgTYkP44wIEZEXYhghj2USJTz3tyK8/00pAGDhmP74wx19W7qqKhUCRvYJl7NEIiJyAoYR8kgNRhPm7ziK3cfKIQjAS5MS8eCIOLnLIiIiF2AYIfmJJqDkIFBbAYREojYqBbO3HcX+n6oQpFRgzdShuOdmrgMhIvJVDCMkL202kLMIMJS1HKpXdEOX+unoHDQSb81Ixq/6dZOxQCIicjWGEZKPNhvYmQG02mkmzFSFrKA1KL1rIG5kECEi8nm8tZfkIZrMMyI2trxTCIAAATd++6J5HBER+TSGEZJHyUGLSzOtCZAAw1nzOCIi8mkMIySP2grnjiMiIq/FMELyCIl07jgiIvJaDCMki63lPVAmhUG0XjJymQCoewBxqe4si4iIZMAwQm4lSRJe+9eP+NPfjmNZUwYEAZDa2vIufQWgULq9RiIici+GEXIbUZSwbJcWq3NPAAAG3PEgMHkLBHWrhmbqGGDyFiBhogxVEhGRu7HPCLlFo1HE0x99h78dNd9B88KEBMwcFQ+gPzDw1xYdWBGXyhkRIiI/wjBCLlfXaMScrUew98Q5BCgErJo8BPcO7XFlgEIJxKfJVyAREcmKYYRc6pe6Rjyy+VscKf0FqkAF1k9Pwu0DIuQui4iIPAjDCDmNSZSQX1yNypp6RISq0CusMx7enI8TFbXQdArExpm3ICnuBrnLJCIiD8MwQk6RU6TDsl1a6PT1LccUAiBKQKQ6GO/NGoH+kaEyVkhERJ6KYYSuW06RDnO2HrHaZaa5h8gTd/ZjECEiojbx1l66LqbLt+u22bsMwBtf/gRT293NiIjIzzGM0HXJL662uDRji05fj/ziajdVRERE3oZhhK5LZU37QcTRcURE5H8YRui6RISqnDqOiIj8DxewUodJkoRvis+3O0YAEKVRISU+zD1FERGR1+nQzMi6desQHx8PlUqFpKQk5OXltTn2k08+wejRo9G9e3eo1WqMHDkSX3zxRYcLJs8gihJe+vtxrPnnjy3H2tjuDksnJECpaP1VIiIiM4fDyI4dOzBv3jwsWbIEhYWFSEtLw7hx41BaWmpz/L59+zB69Gjs3r0bBQUFuOOOOzBhwgQUFhZed/EkD6NJxDMf/wcb9hcDMIeN9dOHI0pjeSkmSqNC1vThSE+MtvU0REREAABBkiSH7rkcMWIEhg8fjqysrJZjgwYNwqRJk5CZmWnXc9x0002YMmUKnn/+ebvGGwwGaDQa6PV6qNVqR8olJ2swmvDU9qPI+b4cSoWAl+8bjPuSegKw7sCaEh/GGREiIj9m7/u3Q2tGGhsbUVBQgGeffdbi+JgxY3Dw4EG7nkMURdTU1CAsrO01BA0NDWhoaGj53GAwOFImucjFBiNmby1A3o9VCFIqsHbaMIy9Karl60qFgJF9wmWskIiIvJFDl2mqqqpgMpkQGRlpcTwyMhLl5eV2PceqVatw8eJFTJ48uc0xmZmZ0Gg0LR+xsbGOlEkuoK9rwvQN3yDvxyp0DlJi08O3WAQRIiKijurQAlZBsJx6lyTJ6pgt27dvxwsvvIAdO3YgIqLtnVsXL14MvV7f8nH69OmOlEkdJZqA4jzg2EdAcR4q9Rcx5a1DKCz9BZpOgdj2uxEY1beb3FUSEZGPcOgyTbdu3aBUKq1mQSorK61mS1rbsWMHZs2ahQ8//BB33313u2ODg4MRHBzsSGnkLNpsIGcRYCi7ckwIR1zDDFSHpuG9WSMwIIr7zBARkfM4NDMSFBSEpKQk5ObmWhzPzc1Fampqm4/bvn07Zs6ciffffx/jx4/vWKXketpsYGeGZRAB0E08j/VBa/D53RcYRIiIyOkcvkyzYMECvPPOO9i4cSOOHz+O+fPno7S0FLNnzwZgvsSSkZHRMn779u3IyMjAqlWrcOutt6K8vBzl5eXQ6/XO+yno+okm84yIjS3vzDfECIg4sNQ8joiIyIkcDiNTpkzBmjVrsHz5cgwdOhT79u3D7t27ERcXBwDQ6XQWPUfefPNNGI1G/OEPf0B0dHTLx1NPPeW8n4KuX8lBqxmRqwmQAMNZ8zgiIiIncrjPiBzYZ8QNjn0EfDzr2uPu2wDcfL/r6yEiIq9n7/s3N8ojs5D2FyA7PI6IiMhODCMEAHhPF4MyKQxim/NkAqDuAcS1vVCZiIioIxhG/JwkSXj93z/iuez/YllTBgQBkNra8i59BaBQur1GIiLybQwjfkySJPx593H8dc8JAMCAOx4EJm+BoG61sZ06Bpi8BUiYKEOVRETk6xxqeka+wyRK+L9PjmHHYXN32+d+nYBZv4oH0B8Y+GvzXTO1FeY1InGpnBEhIiKXYRjxQw1GE+bvOIrdx8qhEIAV9w3G5OSr9v9RKIH4NPkKJCIiv8Iw4mfqGo147L0rO+++9sBQpCdGX/uBRERELsIw4sNMooT84mpU1tQjIlSFAZGh+P17h1FQcgGdApV4KyMJaf26y10mERH5OYYRH5VTpMOyXVro9PUtxwIUAoyiBLUqAJseTkFS3A0yVkhERGTGMOKDcop0mLP1iNUuM8bLTUSevKsfgwgREXkM3trrY0yihGW7tDa2u7tiw/5imNrubkZERORWDCM+Jr+42uLSjC06fT3yi6vdVBEREVH7GEZ8TGVN+0HE0XFERESuxjDiYyJCVU4dR0RE5GpcwOpjLlxsaPfrAoAojQop8WHuKYiIiOgaODPiQz4uOIPHtxe2fN7GdndYOiEBSkXrrxIREcmDYcRHbDn0M/744XcQJeB/k3pi3bThiNJYXoqJ0qiQNX04O64SEZFH4WUaH/DGlz/hL1/8AACYmXojnv91AhQKAWMToyw6sKbEh3FGhIiIPA7DiBeTJAkvf/EDsr46CQB48s6+mD+6PwTBHDiUCgEj+4TLWSIREdE1MYx4KVGUsDT7e7z3dQkAYPG4gXjstj4yV0VEROQ4hhFvIJqAkoNAbQUQEgljz1vxzCff45PCsxAE4KVJiXhwRJzcVRIREXUIw4in02YDOYsAQ1nLIYOyOy5eehBKxQisnjwE9w7tIWOBRERE14dhxJNps4GdGUCrnWa6Gs8hK3ANjo1aiyEMIkRE5OV4a6+nEk3mGREbW94pBEAQBAwpWmEeR0RE5MUYRjxVyUGLSzOtCZAAw1nzOCIiIi/GMOKpaiucO46IiMhDMYx4qpBI544jIiLyUAwjHupH1c2oQDhE6yUjlwmAugcQl+rOsoiIiJyOYcQDHTujx+S38/F84wwIAiC1teVd+gpAoXR7fURERM7EMOJh8ourMe3tr3GhrgnlMaNRN2kTBHWrje3UMcDkLUDCRHmKJCIiciL2GfEge0+cw2PvHUZ9k4gR8WHYMPMWdAkOAAbfa9GBFXGpnBEhIiKfwTDiIXKKdHhieyGaTBLuGNAdWdOToAq8HDgUSiA+Td4CiYiIXIRhxAN8XHAGT3/0HUQJGH9zNF6ZMhRBAbyCRkRE/oFhRGZbDv2M5//2PQBgcnJPZP52MJSK1gtWiYiIfBfDiBuZRAn5xdWorKlHRKgKh0uqsWrPCQDAI6Pi8afxg6BgECEiIj/DMOImOUU6LNulhU5fb/W1p+7qh3l394MgMIgQEZH/YRhxg5wiHeZsPWJjyzuzQdGhDCJEROS3uErSxUyihGW7tG0GEQHAsl1amNputUpEROTTGEZcLL+42ualmWYSAJ2+HvnF1e4rioiIyIMwjLhYZU3bQaQj44iIiHwNw4iLhQbbtywnIlTl4kqIiIg8ExewutCFi4145Z8n2h0jAIjSqJASH+aeooiIiDwMZ0ZcpLKmHlPf+hrHzhoQcnl2pI29d7F0QgIbnRERkd9iGHGBMxfqMHn9IfxQUYNIdTA++0Mq1k8fjiiN5aWYKI0KWdOHIz0xuo1nIiIi8n28TONkxVUX8eDbX6NMX4+eN3TC+7+7Fb3CO6NvRChGJ0RZdGBNiQ/jjAgREfk9hhEn+m+5AdPfyUdVbQP6dO+Cbb+71WI2RKkQMLJPuIwVEhEReR6GESf57vQveGhTPn6pa0JCtBpbZqWgW0iw3GURERF5PIYRR4kmoOQgUFsBhEQCcan45udfMOvdw6htMGJYr67YPDMFms6BcldKRETkFRhGHKHNBnIWAYaylkP1naOwtWYaapuSMbJ3ON55KBld7OwtQkRERLybxn7abGBnhkUQAYCgi+V4VbEaC2N/wKaHb2EQISIichDDiD1Ek3lGxMZ2dwoBEATgDw3vQKV0f2lERETejmHEHiUHrWZEriYAEAxnzeOIiIjIIQwj9qitcO44IiIiasEwYo+QSOeOIyIiohYMI3aQeo2EPjACovWSkcsEQN0DiEt1Z1lEREQ+gWHkGkyihMWfafHMxWkAAKmt7e7SVwAKrmAlIiJyFMNIO5pMIubvOIoPvj2NXCkFXye/AkHdalM7dQwweQuQMFGeIomIiLwcm2K0ob7JhCe2FyJXW4EAhYBXpw5D6uBoQHzIqgMrZ0SIiIg6jmHEhrpGIx7dUoD9P1UhOECB9dOTcMfACPMXFUogPk3eAomIiHwIw0grhvomPLLpWxwuuYAuQUq889At3GmXiIjIhfw2jJhECfnF1aisqUdEqAop8WHQX2pCxsZvUHTWALUqAO8+koJhvW6Qu1QiIiKf1qEFrOvWrUN8fDxUKhWSkpKQl5fX7vi9e/ciKSkJKpUKvXv3xvr16ztUrLPkFOnwq5X/xgNvf42nPjiKB97+GiMz/4Xxr+ah6KwB4V2C8MGjIxlEiIiI3MDhMLJjxw7MmzcPS5YsQWFhIdLS0jBu3DiUlpbaHF9cXIx77rkHaWlpKCwsxP/93//hySefxMcff3zdxXdETpEOc7YegU5fb3G8sqYBOkM9unYOxM7ZI5EQo5alPiIiIn8jSJLUZisvW0aMGIHhw4cjKyur5digQYMwadIkZGZmWo1ftGgRsrOzcfz48ZZjs2fPxnfffYdDhw7Z9T0NBgM0Gg30ej3U6o6HBJMo4Vcr/20VRK4WERqMQ4vvglLRup8IEREROcLe92+HZkYaGxtRUFCAMWPGWBwfM2YMDh60vUncoUOHrMaPHTsWhw8fRlNTk83HNDQ0wGAwWHw4Q35xdbtBBDDPkOQXVzvl+xEREdG1ORRGqqqqYDKZEBlpuQdLZGQkysvLbT6mvLzc5nij0Yiqqiqbj8nMzIRGo2n5iI2NdaTMNlXWtB9EHB1HRERE169DC1gFwfIShiRJVseuNd7W8WaLFy+GXq9v+Th9+nRHyrQSEapy6jgiIiK6fg7d2tutWzcolUqrWZDKykqr2Y9mUVFRNscHBAQgPNx2/47g4GAEBwc7UppdUuLDEK1RoVxfD1sLZQQAURrzbb5ERETkHg7NjAQFBSEpKQm5ubkWx3Nzc5GaanvH2pEjR1qN37NnD5KTkxEYGOhguddHqRCwdEICALS13R2WTkjg4lUiIiI3cvgyzYIFC/DOO+9g48aNOH78OObPn4/S0lLMnj0bgPkSS0ZGRsv42bNno6SkBAsWLMDx48exceNGbNiwAQsXLnTeT+GA9MRoZE0fjiiN5aWYKI0KWdOHIz0xuo1HEhERkSs43IF1ypQpOH/+PJYvXw6dTofExETs3r0bcXFxAACdTmfRcyQ+Ph67d+/G/Pnz8cYbbyAmJgavvfYa7rvvPuf9FA5KT4zG6IQoqw6snBEhIiJyP4f7jMjBWX1GiIiIyH1c0meEiIiIyNkYRoiIiEhWDCNEREQkK4YRIiIikhXDCBEREcmKYYSIiIhkxTBCREREsmIYISIiIlkxjBAREZGsHG4HL4fmJrEGg0HmSoiIiMheze/b12r27hVhpKamBgAQGxsrcyVERETkqJqaGmg0mja/7hV704iiiLKyMoSGhkIQnLeZncFgQGxsLE6fPs09b+zEc+YYni/H8Zw5jufMMTxfjuvoOZMkCTU1NYiJiYFC0fbKEK+YGVEoFOjZs6fLnl+tVvMfpIN4zhzD8+U4njPH8Zw5hufLcR05Z+3NiDTjAlYiIiKSFcMIERERycqvw0hwcDCWLl2K4OBguUvxGjxnjuH5chzPmeN4zhzD8+U4V58zr1jASkRERL7Lr2dGiIiISH4MI0RERCQrhhEiIiKSFcMIERERycrnw8i6desQHx8PlUqFpKQk5OXltTt+7969SEpKgkqlQu/evbF+/Xo3Veo5HDlnn3zyCUaPHo3u3btDrVZj5MiR+OKLL9xYrfwc/TfW7MCBAwgICMDQoUNdW6AHcvScNTQ0YMmSJYiLi0NwcDD69OmDjRs3uqla+Tl6vrZt24YhQ4agc+fOiI6OxsMPP4zz58+7qVr57du3DxMmTEBMTAwEQcBnn312zcf482u/o+fLJa/7kg/74IMPpMDAQOntt9+WtFqt9NRTT0ldunSRSkpKbI4/deqU1LlzZ+mpp56StFqt9Pbbb0uBgYHSRx995ObK5ePoOXvqqaeklStXSvn5+dKJEyekxYsXS4GBgdKRI0fcXLk8HD1fzX755Repd+/e0pgxY6QhQ4a4p1gP0ZFzNnHiRGnEiBFSbm6uVFxcLH3zzTfSgQMH3Fi1fBw9X3l5eZJCoZBeffVV6dSpU1JeXp500003SZMmTXJz5fLZvXu3tGTJEunjjz+WAEiffvppu+P9/bXf0fPlitd9nw4jKSkp0uzZsy2ODRw4UHr22Wdtjn/mmWekgQMHWhx77LHHpFtvvdVlNXoaR8+ZLQkJCdKyZcucXZpH6uj5mjJlivSnP/1JWrp0qd+FEUfP2T/+8Q9Jo9FI58+fd0d5HsfR8/WXv/xF6t27t8Wx1157TerZs6fLavRk9ry58rX/CnvOly3X+7rvs5dpGhsbUVBQgDFjxlgcHzNmDA4ePGjzMYcOHbIaP3bsWBw+fBhNTU0uq9VTdOSctSaKImpqahAWFuaKEj1KR8/Xpk2bcPLkSSxdutTVJXqcjpyz7OxsJCcn4+WXX0aPHj3Qv39/LFy4EJcuXXJHybLqyPlKTU3FmTNnsHv3bkiShIqKCnz00UcYP368O0r2Sv7+2n+9nPG67xUb5XVEVVUVTCYTIiMjLY5HRkaivLzc5mPKy8ttjjcajaiqqkJ0dLTL6vUEHTlnra1atQoXL17E5MmTXVGiR+nI+frxxx/x7LPPIi8vDwEBPvvfr00dOWenTp3C/v37oVKp8Omnn6Kqqgpz585FdXW1z68b6cj5Sk1NxbZt2zBlyhTU19fDaDRi4sSJWLt2rTtK9kr+/tp/vZzxuu+zMyPNBEGw+FySJKtj1xpv67gvc/ScNdu+fTteeOEF7NixAxEREa4qz+PYe75MJhOmTZuGZcuWoX///u4qzyM58m9MFEUIgoBt27YhJSUF99xzD1avXo3Nmzf7xewI4Nj50mq1ePLJJ/H888+joKAAOTk5KC4uxuzZs91Rqtfia3/HOOt132d/NevWrRuUSqXVbw+VlZVWCbhZVFSUzfEBAQEIDw93Wa2eoiPnrNmOHTswa9YsfPjhh7j77rtdWabHcPR81dTU4PDhwygsLMTjjz8OwPxGK0kSAgICsGfPHtx5551uqV0uHfk3Fh0djR49elhsQz5o0CBIkoQzZ86gX79+Lq1ZTh05X5mZmRg1ahSefvppAMDgwYPRpUsXpKWl4aWXXuJv+Tb4+2t/Rznzdd9nZ0aCgoKQlJSE3Nxci+O5ublITU21+ZiRI0dajd+zZw+Sk5MRGBjoslo9RUfOGWBOxjNnzsT777/vV9elHT1farUax44dw9GjR1s+Zs+ejQEDBuDo0aMYMWKEu0qXTUf+jY0aNQplZWWora1tOXbixAkoFAr07NnTpfXKrSPnq66uDgqF5Uu7UqkEcOW3fbLk76/9HeH01/0OL331As23xG3YsEHSarXSvHnzpC5dukg///yzJEmS9Oyzz0ozZsxoGd98e9f8+fMlrVYrbdiwwa9u75Ikx8/Z+++/LwUEBEhvvPGGpNPpWj5++eUXuX4Et3L0fLXmj3fTOHrOampqpJ49e0r333+/9P3330t79+6V+vXrJ/3ud7+T60dwK0fP16ZNm6SAgABp3bp10smTJ6X9+/dLycnJUkpKilw/gtvV1NRIhYWFUmFhoQRAWr16tVRYWNhyOzRf+y05er5c8brv02FEkiTpjTfekOLi4qSgoCBp+PDh0t69e1u+9tBDD0m33XabxfivvvpKGjZsmBQUFCTdeOONUlZWlpsrlp8j5+y2226TAFh9PPTQQ+4vXCaO/hu7mj+GEUly/JwdP35cuvvuu6VOnTpJPXv2lBYsWCDV1dW5uWr5OHq+XnvtNSkhIUHq1KmTFB0dLT344IPSmTNn3Fy1fL788st2X5f42m/J0fPlitd9QZI4b0dERETy8dk1I0REROQdGEaIiIhIVgwjREREJCuGESIiIpIVwwgRERHJimGEiIiIZMUwQkRERLJiGCEiIiJZMYwQERGRrBhGiIiISFYMI0RERCQrhhEiIiKS1f8Hkzb2XEAyUqAAAAAASUVORK5CYII=", + "text/plain": [ + "
" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "plt.plot(domain_x, f(domain_x), '-o')\n", + "plt.plot(x_, f_x, 'o')\n", + "plt.legend(['original', 'half'])" + ] + }, + { + "cell_type": "code", + "execution_count": 11, + "id": "0c6a644b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Riemann sum integral: 0.6170376421171327\n", + "Exact Integral in first domain: 0.6173165676349102\n" + ] + } + ], + "source": [ + "Riemann = (np.sum(f_x)*(b-a))/2**n\n", + "print(\"Riemann sum integral: {}\".format(Riemann))\n", + "print('Exact Integral in first domain: {}'.format(sin_integral(a,b)))" + ] + }, + { + "cell_type": "markdown", + "id": "785e8e08", + "metadata": {}, + "source": [ + "### 2.3. Array Normalisation\n", + "\n", + "The idea is to encode the $2^n$ discretized array $f_{x_i}$ in a $n+1$ qubit circuit. Before doing that we are going to normalise the array in the following way:\n", + "\n", + "\n", + "$$f^{norm}_{x_i} = \\frac{f_{x_i}}{\\max(|f_{x_i}|)}$$\n", + "\n", + "If $\\max{|f_{x_i}|} \\leq 1$ then this step can be omitted.\n", + "\n", + "Now the computed integral will be\n", + "\n", + "\\begin{equation}\n", + "S_{[a,b]} = \\frac{b-a}{2^n} \\sum_{i=0}^{2^n-1} f_{x_i} = \\frac{b-a}{2^n} \\sum_{i=0}^{2^n-1} \\max(|f_{x_i}|) f^{norm}_{x_i} =\\frac{\\max(|f_{x_i}|)(b-a)}{2^n} \\sum_{i=0}^{2^n-1} f^{norm}_{x_i}\n", + "\\end{equation}" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "092a27f7", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Normalization constant: 0.9085519168534011\n", + "Riemman sum integral: 0.6170376421171329\n", + "Exact Integral in first domain: 0.6173165676349102\n" + ] + } + ], + "source": [ + "normalization = np.max(np.abs(f_x))\n", + "print(\"Normalization constant: {}\".format(normalization))\n", + "#normalization = 1.0\n", + "f_norm_x = f_x/normalization\n", + "Riemann = normalization * np.sum(f_norm_x)*(b-a)/2**n\n", + "#Now we need to be aware of the normalization constant when computing Rieman sum\n", + "print(\"Riemman sum integral: {}\".format(Riemann))\n", + "print('Exact Integral in first domain: {}'.format(sin_integral(a,b)))" + ] + }, + { + "cell_type": "markdown", + "id": "3ca9da3d-e350-43a4-ac33-8aa6c3fd2ebf", + "metadata": {}, + "source": [ + "### 2.4. Encoding function in a quantum circuit.\n", + "\n", + "The next step is to codify $f^{norm}_{x_i}$ array in a quantum circuit. The following procedure must be used:\n", + "\n", + "1. Initialize a quantum register with at least $n+1$ qubits, where $n$ must be equal to the $n$ used to define the $2^n$ discretization intervals $$|0\\rangle\\otimes|0\\rangle_n \\tag{1}$$\n", + "2. Apply the uniform distribution over the first $n$ qubits: $$\\big(I \\otimes H^{\\otimes n}\\big)\\big(|0\\rangle \\otimes|0\\rangle_{n}\\big) = |0\\rangle \\otimes H^{\\otimes n}|0\\rangle_{n}=\n", + "\\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1}|0\\rangle \\otimes|i\\rangle_{n} \\tag{2}$$\n", + "3. Creates an operator $\\mathbf{U}_f$ for encoding the $f\\_norm_{x_i}$. This operator acts in the following way: $$\\mathbf{U}_f \\left(|0\\rangle \\otimes |i\\rangle_n\\right) = \\big(f^{norm}_{x_i}|0\\rangle + \\beta_i |1\\rangle \\big) \\otimes |i\\rangle_n \\tag{3}$$ In the equation $(3)$ the coefficient of $|1\\rangle$ it is not important for us the only important coefficient it is the $|0\\rangle$ one.\n", + "4. Apply the $\\mathbf{U}_f$ operator over $n+1$ qubits: $$\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n}\\right)|0\\rangle\\otimes|0\\rangle_{n} \\tag{4}$$\n", + "5. Applying equation $(2)$ and $(3)$ on $(4)$: $$\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\mathbf{U}_f \\left(\\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} |0\\rangle\\otimes|i\\rangle_{n}\\right) = \\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} \\mathbf{U}_f \\left(|0\\rangle\\otimes|i\\rangle_{n}\\right) = \\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} |i\\rangle_{n} \\otimes \\left(f^{norm}_{x_i}|0\\rangle + \\beta_i|1\\rangle \\right) \\tag{5}$$\n", + "6. Finally the uniform distribution will be applied over the first n qubits again: $$|\\Psi \\rangle = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} \\tag{6}$$\n", + "7. So applying $(5)$ on $(6)$: $$|\\Psi \\rangle = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} H^{\\otimes n}|i\\rangle_{n} \\otimes \\left(f^{norm}_{x_i}|0\\rangle + \\beta_i|1\\rangle \\right) \\tag{7}$$\n", + "8. We are interested only in $|0\\rangle \\otimes |i\\rangle_{n}$ so we don't need to take into account other terms, so $(7)$ can be expressed as: $$|\\Psi \\rangle = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\frac{1}{\\sqrt{2^n}} \\sum_{i=0}^{2^{n}-1} f^{norm}_{x_i}|0\\rangle \\otimes H^{\\otimes n}|i\\rangle_{n} + \\cdots \\tag{8}$$\n", + "9. It is known that: $$H^{\\otimes n} = \\frac{1}{\\sqrt{2^n}} \\sum_{j=0}^{2^n}\\sum_{k=0}^{2^n} (-1)^{jk} |j\\rangle_n {}_{n} \\langle k|$$.\n", + "10. So:$$H^{\\otimes n} |i\\rangle_n = \\frac{1}{\\sqrt{2^n}} \\sum_{j=0}^{2^n}\\sum_{k=0}^{2^n} (-1)^{jk} |j\\rangle_n {}_{n} \\langle k|i\\rangle_n = \\frac{1}{\\sqrt{2^n}} \\sum_{j=0}^{2^n} (-1)^{ji} |j\\rangle = \\frac{1}{\\sqrt{2^n}} |0\\rangle_n + \\frac{1}{\\sqrt{2^n}} \\sum_{j=1}^{2^n} (-1)^{ji} |j\\rangle \\tag{9}$$\n", + "11. Finally applying $(9)$ in $(8)$ and taking only into account the $|0\\rangle \\otimes |0\\rangle_{n}$ $$|\\Psi \\rangle = \\left(I\\otimes H^{\\otimes n} \\right)\\mathbf{U}_f\\left(I\\otimes H^{\\otimes n} \\right)|0\\rangle\\otimes|0\\rangle_{n} = \\frac{1}{2^n} \\sum_{i=0}^{2^{n}-1} f^{norm}_{x_i} |0\\rangle\\otimes|0\\rangle_{n} + \\cdots \\tag{10}$$\n", + "\n", + "Using this procedure the two mandatory operators $\\mathbf{A}^I(f_{x_i})$ (one for each integration interval) can be created:\n", + "$$\\mathbf{A}^I(f_{x_i}) = \\left(\\mathbb{I}\\otimes H^{\\otimes n} \\right)\\mathbf{U}^I_f\\left(\\mathbb{I}\\otimes H^{\\otimes n} \\right)$$\n", + "\n", + "Using $(10)$ the behaviour of these operators can be summarised as:\n", + "\n", + "$$|\\Psi \\rangle = \\mathbf{A}^I(f_{x_i}) |0\\rangle\\otimes|0\\rangle_{n}= \\frac{1}{2^n} \\sum_{i=0}^{2^{n}-1} f^{norm}_{x_i} |0\\rangle\\otimes|0\\rangle_{n} + \\cdots \\tag{11}$$\n", + "\n", + "$(11)$ can be compared with the equation of the **AE kernel** (see notebook *01_AmplitudeEstimationKernel.ipynb*):\n", + "\n", + "$$|\\Psi\\rangle= \\mathbf{A}|0\\rangle_n = \\sqrt{a} |\\Psi_0\\rangle + \\sqrt{1-a}|\\Psi_1\\rangle$$\n", + "\n", + "Where $|\\Psi_0\\rangle = |0\\rangle\\otimes|0\\rangle_{n}$ and $$\\sqrt{a} = \\frac{1}{2^n} \\sum_{i=0}^{2^{n}-1} f^{norm}_{x_i}$$\n", + "\n", + "Now the Riemann sum approximation of the desired integral can be computed by measuring the probability of obtaining the state $|\\Psi_0\\rangle = |0\\rangle \\otimes |0\\rangle_n$ as shown in $(12)$:\n", + "$$\\mathbf{P}[|\\Psi_0\\rangle] = \\left| \\; \\langle \\Psi_0\\ |\\Psi\\rangle \\; \\right|^2 = \\left| \\; \\langle \\Psi_0\\ | \\frac{1}{2^n} \\sum_{i=0}^{2^{n}-1} f^{norm}_{x_i} |\\Psi_0\\rangle\\; \\right|^2 = \\left| \\frac{1}{2^n} \\sum_{i=0}^{2^{n}-1} f^{norm}_{x_i} \\right|^2 = \\tilde{a} \\tag{12}$$\n", + "\n", + "The $\\thicksim$ in $\\tilde{a}$ indicates that the amplitude was obtained using a quantum measurement. Now plugging $(12)$ into the definition of the integral:\n", + "\n", + "$$\\tilde{S}_{[a,b]}= \\frac{\\max(f_{x_i}) \\left(b-a\\right)}{2^n}\\left(2^n\\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}\\right) \\tag{13}$$\n", + "\n", + "The $\\thicksim$ in $\\tilde{S}_{[a,b]}$ indicates that the integral was obtained using a measurement meanwhile the $S_{[a,b]}$ is for pure Riemann sum calculation as shown" + ] + }, + { + "cell_type": "markdown", + "id": "bf0e3a5d", + "metadata": {}, + "source": [ + "#### Operator $\\mathbf{U}_f$\n", + "\n", + "Here we are going to explain how to build the operator $\\mathbf{U}_f$:\n", + "\n", + "1. Given the array $f^{norm}_{x_i}$ we need to compute: $\\phi_{x_i} = \\arccos({f^{norm}_{x_i}})$.\n", + "2. for a given state $|i\\rangle_n \\otimes |0\\rangle$ we need to implement a rotation around the *y-axis* over the last qubit ($|0\\rangle$) controlled by the state $|i\\rangle_n$ of $2*\\phi_{x_i}$. So we need to build following operation:\n", + "\n", + "$$|i\\rangle_n \\otimes |0\\rangle \\rightarrow |i\\rangle_n \\otimes \\mathbf{R}_y(2*\\phi_{x_i})|0\\rangle = |i\\rangle_n \\otimes \\left( \\cos(\\phi_{x_i})|0\\rangle + \\sin(\\phi_{x_i})|1\\rangle \\right) $$\n", + "\n", + "3. Now undoing the $\\phi_{x_i}$ and doing $\\beta_i = \\sin(\\phi_{x_i})$ we can obtain the desired operator $\\mathbf{U}_f$:\n", + "\n", + "$$|i\\rangle_n \\otimes \\left(f^{norm}_{x_i} |0\\rangle + \\beta_i|1\\rangle \\right) = \\mathbf{U}_f |i\\rangle_n \\otimes |0\\rangle$$" + ] + }, + { + "cell_type": "markdown", + "id": "ca85032e", + "metadata": {}, + "source": [ + "The encoding of the array $f^{norm}_{x_i}$ in a unitary operator $\\mathbf{A}^I(f_{x_i})$ is done transparently by python class *Encoding* from the *QQuantLib/DL/encoding_protocols* module.\n", + "\n", + "This class creates the operator $\\mathbf{A}^I(f_{x_i})$ under the *oracle* property of the *Encoding* class. \n", + "\n", + "For instantiating the class the following arguments should be provided:\n", + "\n", + "* array_function: the array function $f^{norm}_{x_i}$.\n", + "* array_probability: should be fixed to None\n", + "* encoding: should be fixed to 2 (the class implement until 3 different encoding protocols but for the **BTC** the it should be fixed to 2)" + ] + }, + { + "cell_type": "code", + "execution_count": 15, + "id": "a01a0395", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.path.append(\"../../\")" + ] + }, + { + "cell_type": "code", + "execution_count": 16, + "id": "193f009a", + "metadata": {}, + "outputs": [], + "source": [ + "from QQuantLib.DL.encoding_protocols import Encoding" + ] + }, + { + "cell_type": "code", + "execution_count": 17, + "id": "49386077", + "metadata": {}, + "outputs": [], + "source": [ + "encoding_object = Encoding(\n", + " array_function=f_norm_x, \n", + " array_probability=None, \n", + " encoding=2\n", + ")\n", + "encoding_object.run()" + ] + }, + { + "cell_type": "code", + "execution_count": 18, + "id": "03299414", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "Q0Q1Q2Q3Q4UD [4]F_{Function}UD [4]†" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 18, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "operator_A = encoding_object.oracle\n", + "%qatdisplay operator_A --svg" + ] + }, + { + "cell_type": "markdown", + "id": "440db60f", + "metadata": {}, + "source": [ + "Additionally the *function_gate* attribute give us acces to the operator: $\\mathbf{U}_g$" + ] + }, + { + "cell_type": "code", + "execution_count": 19, + "id": "a08cc3b0", + "metadata": {}, + "outputs": [ + { + "data": { + "text/html": [ + "Q0Q1Q2Q3Q4RY [ 1.76]RY [ 0.73]RY [- 0.06]RY [ 0.39]RY [- 0.04]RY [ 0.03]RY [- 0.04]RY [ 0.21]RY [- 0.02]RY [ 0.02]RY [- 0.02]RY [ 0.02]RY [- 0.03]RY [ 0.02]RY [- 0.03]RY [ 0.11]" + ], + "text/plain": [ + "" + ] + }, + "execution_count": 19, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "operator_Uf = encoding_object.function_gate\n", + "%qatdisplay operator_Uf --depth --svg" + ] + }, + { + "cell_type": "markdown", + "id": "36788404", + "metadata": {}, + "source": [ + "## 7. Grover-like operator of $\\mathbf{A}$\n", + "\n", + "For each one of the 3 $\\mathbf{A(f_{x_i})}$ operators the Grover-like operator should be constructed using:\n", + "\n", + "$$\\mathbf{G}(\\mathbf{A(f_{x_i})}) = \\mathbf{A(f_{x_i})} \\left(\\hat{I} - 2|0\\rangle\\langle 0|\\right) \\mathbf{A(f_{x_i})}^{\\dagger}\\left(\\hat{I} - 2|\\Psi_0\\rangle\\langle \\Psi_0|\\right)$$" + ] + }, + { + "cell_type": "markdown", + "id": "94318e62", + "metadata": {}, + "source": [ + "In the case of our **QQuantLib** the Grover-like operator building can be done in an easy way using the **grover** function from *QQuantLib.AA.amplitude_amplification* module. The input of this function will be the following attribtues of the encoding class:\n", + "\n", + "* *oracle*\n", + "* *target*\n", + "* *index*" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f42422c5", + "metadata": {}, + "outputs": [], + "source": [ + "from QQuantLib.AA.amplitude_amplification import grover" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b6d320c1", + "metadata": {}, + "outputs": [], + "source": [ + "operator_G = grover(\n", + " oracle = encoding_object.oracle,\n", + " target = encoding_object.target,\n", + " index = encoding_object.index\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bac5d25a", + "metadata": {}, + "outputs": [], + "source": [ + "%qatdisplay operator_G --depth 2 --svg" + ] + }, + { + "cell_type": "markdown", + "id": "46e09369", + "metadata": {}, + "source": [ + "## 8. Amplitude Estimation Algorithm\n", + "\n", + "For each of the 3 $\\mathbf{A(f_{x_i})}$ operators and their correspondent Grover-like operators $\\mathbf{G}(\\mathbf{A(f_{x_i})})$ the desired *Amplitude Estimation* algorithm can be used. \n", + "\n", + "An amplitude estimation algorithm, usually, returns the probabiliy of getting the state $|\\Psi_0\\rangle$, so a post post-proccesing, for getting the estimator of the integral, should be done, as explained in section 6:\n", + "\n", + "$$\\tilde{S}_{[a,b]} = \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}\\right)$$\n", + "\n", + "\n", + "In case the algorithm returns some different value an additional post-processing should be done in order to get the $\\tilde{S}_{[a,b]}$ integral estimator." + ] + }, + { + "cell_type": "markdown", + "id": "794eee23", + "metadata": {}, + "source": [ + "In the case of **QQuantLib** steps 7 and 8 can be done straightforward using the different *Amplitude Estimation* classes of the modules from package *QQuantLib.AE*. \n", + "In fact using the **AE** class from *QQuantLib.AE.ae_class* we can acces to all the implemented algorithms in an easy way" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e84ed34a", + "metadata": {}, + "outputs": [], + "source": [ + "#This cell loads the QLM solver.\n", + "#QLMaaS == False -> uses PyLinalg\n", + "#QLMaaS == True -> try to use LinAlg (for using QPU as CESGA QLM one)\n", + "from QQuantLib.utils.qlm_solver import get_qpu\n", + "linalg_qpu = get_qpu('c')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6db22fa8", + "metadata": {}, + "outputs": [], + "source": [ + "from QQuantLib.AE.ae_class import AE" + ] + }, + { + "cell_type": "markdown", + "id": "e16b209a", + "metadata": {}, + "source": [ + "First we create a base configuration dictionary for the **AE** object" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "57a5da80", + "metadata": {}, + "outputs": [], + "source": [ + "#AE base configuration dictionary\n", + "ae_dictionary = {\n", + " #QPU\n", + " 'qpu': linalg_qpu,\n", + " #Multi controlled decomposition\n", + " 'mcz_qlm': False, \n", + " \n", + " #shots\n", + " 'shots': None,\n", + " \n", + " #MLAE\n", + " 'schedule': [\n", + " [],\n", + " []\n", + " ],\n", + " 'delta' : None,\n", + " 'ns' : None,\n", + " \n", + " #CQPEAE\n", + " 'auxiliar_qbits_number': None,\n", + " #IQPEAE\n", + " 'cbits_number': None,\n", + " #IQAE & RQAQE\n", + " 'epsilon': None,\n", + " #IQAE\n", + " 'alpha': None,\n", + " #RQAE\n", + " 'gamma': None,\n", + " 'q': None\n", + "}" + ] + }, + { + "cell_type": "markdown", + "id": "349eda8e", + "metadata": {}, + "source": [ + "We are going to use the **IQAE** algorithm for solving the integral" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "73c20e30", + "metadata": {}, + "outputs": [], + "source": [ + "ae_dictionary.update({\n", + " 'ae_type': 'IQAE',\n", + " 'epsilon' : 0.001,\n", + " 'alpha': 0.05\n", + "})\n", + "ae_object = AE(\n", + " oracle=encoding_object.oracle,\n", + " target=encoding_object.target,\n", + " index=encoding_object.index,\n", + " **ae_dictionary\n", + ")\n", + "ae_object.run()\n", + "#We need to post-procces the return in order to get the correct estimator\n", + "ae_estimator_S = (b-a)*normalization*(2**n*np.sqrt(ae_object.ae_pdf))/2**n" + ] + }, + { + "cell_type": "markdown", + "id": "8e61b5e8", + "metadata": {}, + "source": [ + "In the case of the **IQAE** algorithm the lower and te upper limits for the amplitude is provided!!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6e1d389e", + "metadata": {}, + "outputs": [], + "source": [ + "ae_estimator_S" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "21698868", + "metadata": {}, + "outputs": [], + "source": [ + "print('Integral in first domain: {}'.format(sin_integral(a,b)))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "a60e682a", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Is it the integral between the bound limits of the estimator: {}\".format(\n", + " (sin_integral(a,b) < ae_estimator_S['ae_u'])[0] & (sin_integral(a,b) > ae_estimator_S['ae_l'])[0]\n", + "))\n" + ] + }, + { + "cell_type": "markdown", + "id": "197ae2a6", + "metadata": {}, + "source": [ + "In the case of the **RQAE** algorithm the returned value is directly the **AMplitude** and not he probability (like in the ofhter methods implemented in the libary) so an additionall post-procces is mandatory. In this case:\n", + "\n", + "$$\\tilde{S}_{[a,b]} = \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\mathbf{P}[|\\Psi_0\\rangle]\\right)$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "28107426", + "metadata": {}, + "outputs": [], + "source": [ + "ae_dictionary.update({\n", + " 'ae_type': 'RQAE',\n", + " 'epsilon' : 0.001,\n", + " 'alpha': None,\n", + " 'gamma': 0.05\n", + "})\n", + "ae_object = AE(\n", + " oracle=encoding_object.oracle,\n", + " target=encoding_object.target,\n", + " index=encoding_object.index,\n", + " **ae_dictionary\n", + ")\n", + "ae_object.run()\n", + "\n", + "#We need to post-procces the return in order to get the correct estimator.\n", + "#In the RQAE case the amplitude instead of the probability is returned!!!\n", + "ae_estimator_S = (b-a)*normalization*2**n*ae_object.ae_pdf/2**n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b29f2100", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Is it the integral between the bound limits of the estimator: {}\".format(\n", + " (sin_integral(a,b) < ae_estimator_S['ae_u'])[0] & (sin_integral(a,b) > ae_estimator_S['ae_l'])[0]\n", + "))\n" + ] + }, + { + "cell_type": "markdown", + "id": "f3891ada", + "metadata": {}, + "source": [ + "#### q_solve_integral\n", + "\n", + "The *q_solve_integral* function from *QQuantLib.finance.quantum_integration* module allows to solve integrals by *Amplitude Estimation* techniques of inputs arrays. \n", + "\n", + "The input of the function will be a complete dictionary with the configuration of the *amplitude estimation* algorithm and information about the arrays:\n", + "\n", + "* array_function : numpy array with the desired array function for Riemann sum\n", + "* array_probability : numpy array a probability distribtuion for computation of expected values. In our case this will be None.\n", + "* encoding : int for selecting the encoding. In the case of the benchmark will be 2.\n", + "\n", + "The outputs of the *q_solve_integral* function are:\n", + "* ae_estimation: pandas DataFrame with the desired integral computation and the upper and lower limits if applied.\n", + "* solver_ae: objet based on the AE class\n", + "\n", + "The *q_solve_integral* returns directly the integral of the input function directly. So the returned will be: $2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}$ (or the $2^{n} \\mathbf{P}[|\\Psi_0\\rangle]$ for the **RQAE** algorithm)\n", + "\n", + "\n", + "**BE AWARE!!**\n", + "\n", + "This is why we kept the $2^n$ in the integral computation, for taking advance of the implemented integral in the *q_solve_integral* function!!\n", + "\n", + "$$\n", + "\\tilde{S}_{[a,b]} = \\left\\{\n", + "\\begin{array}{ll}\n", + " \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\sqrt{\\mathbf{P}[|\\Psi_0\\rangle]}\\right) & For \\; probability \\; measurements \\\\\n", + " \\frac{\\max(f_{x_i}) (b-a)}{2^n} \\left(2^{n} \\mathbf{P}[|\\Psi_0\\rangle]\\right) & For \\; amplitude \\; measurements \\\\\n", + "\\end{array} \n", + "\\right.\n", + "$$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "868865f1", + "metadata": {}, + "outputs": [], + "source": [ + "from QQuantLib.finance.quantum_integration import q_solve_integral" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ca9d8a30", + "metadata": {}, + "outputs": [], + "source": [ + "#AE base configuration dictionary\n", + "ae_dictionary = {\n", + " #QPU\n", + " 'qpu': linalg_qpu,\n", + " #Multi controlled decomposition\n", + " 'mcz_qlm': False, \n", + " \n", + " #shots\n", + " 'shots': None,\n", + " \n", + " #MLAE\n", + " 'schedule': [\n", + " [],\n", + " []\n", + " ],\n", + " 'delta' : None,\n", + " 'ns' : None,\n", + " \n", + " #CQPEAE\n", + " 'auxiliar_qbits_number': None,\n", + " #IQPEAE\n", + " 'cbits_number': None,\n", + " #IQAE & RQAQE\n", + " 'epsilon': None,\n", + " #IQAE\n", + " 'alpha': None,\n", + " #RQAE\n", + " 'gamma': None,\n", + " 'q': None\n", + "}\n", + "#IQAE configuration\n", + "ae_dictionary.update({\n", + " 'ae_type': 'IQAE',\n", + " 'epsilon' : 0.001,\n", + " 'alpha': 0.05\n", + "})\n", + "\n", + "encoding_dict = {\n", + " \"array_function\" : f_norm_x,\n", + " \"array_probability\" : None,\n", + " \"encoding\" : 2 \n", + "}\n", + "\n", + "ae_dictionary.update(encoding_dict)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "99792f3d", + "metadata": {}, + "outputs": [], + "source": [ + "solution, solver_object = q_solve_integral(**ae_dictionary)\n", + "#Integral of the input array is returned so only normalization has to be took into account\n", + "ae_estimator_S = normalization*solution*(b-a)/2**n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "e3f374ef", + "metadata": {}, + "outputs": [], + "source": [ + "ae_estimator_S" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "bebd04a9", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Is it the integral between the bound limits of the estimator: {}\".format(\n", + " (sin_integral(a,b) < ae_estimator_S['ae_u'])[0] & (sin_integral(a,b) > ae_estimator_S['ae_l'])[0]\n", + "))" + ] + }, + { + "cell_type": "markdown", + "id": "4ea9a2d5", + "metadata": {}, + "source": [ + "When using the **RQAE** algorithm *q_solve_integral* deals with the internal normalisations and the integral of the input array is returned in an transparent way." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cad85df4", + "metadata": {}, + "outputs": [], + "source": [ + "#AE base configuration dictionary\n", + "ae_dictionary = {\n", + " #QPU\n", + " 'qpu': linalg_qpu,\n", + " #Multi controlled decomposition\n", + " 'mcz_qlm': False, \n", + " \n", + " #shots\n", + " 'shots': None,\n", + " \n", + " #MLAE\n", + " 'schedule': [\n", + " [],\n", + " []\n", + " ],\n", + " 'delta' : None,\n", + " 'ns' : None,\n", + " \n", + " #CQPEAE\n", + " 'auxiliar_qbits_number': None,\n", + " #IQPEAE\n", + " 'cbits_number': None,\n", + " #IQAE & RQAQE\n", + " 'epsilon': None,\n", + " #IQAE\n", + " 'alpha': None,\n", + " #RQAE\n", + " 'gamma': None,\n", + " 'q': None\n", + "}\n", + "#IQAE configuration\n", + "ae_dictionary.update({\n", + " 'ae_type': 'RQAE',\n", + " 'epsilon' : 0.001,\n", + " 'gamma': 0.05,\n", + " 'q' : 1.2\n", + "})\n", + "\n", + "encoding_dict = {\n", + " \"array_function\" : f_norm_x,\n", + " \"array_probability\" : None,\n", + " \"encoding\" : 2 \n", + "}\n", + "\n", + "ae_dictionary.update(encoding_dict)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "ac004a01", + "metadata": {}, + "outputs": [], + "source": [ + "solution, solver_object = q_solve_integral(**ae_dictionary)\n", + "#Integral of the input array is returned so only normalization has to be took into account\n", + "ae_estimator_S = normalization*solution*(b-a)/2**n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c778e480", + "metadata": {}, + "outputs": [], + "source": [ + "ae_estimator_S" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "c6be726d", + "metadata": {}, + "outputs": [], + "source": [ + "print(\"Is it the integral between the bound limits of the estimator: {}\".format(\n", + " (sin_integral(a,b) < ae_estimator_S['ae_u'])[0] & (sin_integral(a,b) > ae_estimator_S['ae_l'])[0]\n", + "))" + ] + }, + { + "cell_type": "markdown", + "id": "9c1a0dd0", + "metadata": {}, + "source": [ + "## 9. Getting the metrics\n", + "\n", + "Once the amplitude estimator of the integral, $\\tilde{S}$, is obtained following metrics should be computed:\n", + "\n", + "* *Absolute error* between the ae estimator and the exact integral to compute: $\\epsilon = |\\tilde{S} - (\\cos(a)-\\cos(b))|$\n", + "\n", + "* *Oracle calls*: total number of calls of the operator $\\mathbf{A}$ (the number of shots should be taking into account in this calculation).\n", + "\n", + "* *Elapsed time*:The complete time for getting $\\tilde{S}$. This time will include all the complete steps this is (it is not necesary have times of each individual step):\n", + " * Discretization time\n", + " * Creation of operator $\\mathbf{A}$\n", + " * Creation of the Grover-like operator: $\\mathbf{G}$\n", + " * Execution of the *amplitude Estimation* algorithm\n", + " * Post-processing execution\n", + " * Metrics calculation." + ] + }, + { + "cell_type": "markdown", + "id": "9e7a56cb", + "metadata": {}, + "source": [ + "\n", + "The *sine_integral* function from *AmplitudeEstimation/ae_sine_integral* allows the user perform one complete computation of the sine integral by **AE** techniques. The inputs are:\n", + "\n", + "* *n_qbits* : number of qubits for interval discretization.\n", + "* *interval* : int for selecting the integration interval.\n", + "* *ae_dictionary* : python dictionary with the complete configuration of the *Amplitude Estimation* algorithm.\n", + "\n", + "The return is a panda DataFrames with the complete information of the benchamrk.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3a802720", + "metadata": {}, + "outputs": [], + "source": [ + "#AE base configuration dictionary\n", + "ae_dictionary = {\n", + " #QPU\n", + " 'qpu': 'c',\n", + " #Multi controlled decomposition\n", + " 'mcz_qlm': False, \n", + " \n", + " #shots\n", + " 'shots': None,\n", + " \n", + " #MLAE\n", + " 'schedule': [\n", + " [],\n", + " []\n", + " ],\n", + " 'delta' : None,\n", + " 'ns' : None,\n", + " \n", + " #CQPEAE\n", + " 'auxiliar_qbits_number': None,\n", + " #IQPEAE\n", + " 'cbits_number': None,\n", + " #IQAE & RQAQE\n", + " 'epsilon': None,\n", + " #IQAE\n", + " 'alpha': None,\n", + " #RQAE\n", + " 'gamma': None,\n", + " 'q': None\n", + "}\n", + "#IQAE configuration\n", + "ae_dictionary.update({\n", + " 'ae_type': 'IQAE',\n", + " 'epsilon' : 0.001,\n", + " 'alpha': 0.05\n", + "})" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "2b30ff7e", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "from ae_sine_integral import sine_integral" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "170d2a4b", + "metadata": {}, + "outputs": [], + "source": [ + "#first interval integral benchmarking\n", + "pdf = sine_integral(4, 0, ae_dictionary)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cb94f6da", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "pdf[['integral_ae', 'integral_ae_l', 'integral_ae_u', 'exact_integral', 'absolute_error_sum']]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cf95c732", + "metadata": {}, + "outputs": [], + "source": [ + "#second integral is negative\n", + "pdf = sine_integral(4, 1, ae_dictionary)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "cd7df328", + "metadata": { + "scrolled": true + }, + "outputs": [], + "source": [ + "pdf[['integral_ae', 'integral_ae_l', 'integral_ae_u', 'exact_integral', 'absolute_error_sum']]" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "3afd06ee", + "metadata": {}, + "outputs": [], + "source": [ + "# Only RQAE can deal with negative integrals\n", + "#RQAE configuration\n", + "ae_dictionary.update({\n", + " 'ae_type': 'RQAE',\n", + " 'epsilon' : 0.001,\n", + " 'gamma': 0.05,\n", + " 'q' : 2.0,\n", + " 'alpha': None\n", + "})\n", + "#second integral is negative\n", + "pdf = sine_integral(4, 1, ae_dictionary)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "b9be8c50", + "metadata": {}, + "outputs": [], + "source": [ + "pdf[['integral_ae', 'integral_ae_l', 'integral_ae_u', 'exact_integral', 'absolute_error_sum']]" + ] + }, + { + "cell_type": "markdown", + "id": "d85f36ed", + "metadata": {}, + "source": [ + "## 10. Benchmark Summary\n", + "\n", + "We have all the pieces needed for do a complete benchmarking. Now we summarize the procedure:\n", + "\n", + "1. Execute following steps for each one of the 3 integration described intervals in Sub-Section 6.1\n", + " 1. Execute the following steps from n = 4 to the maximum nubmer of qubits. For each n execute following steps 10 times.\n", + " 1. Create the domain discretization (sub-section 6.2)\n", + " 2. Create the array with the correspondient sine function discretization (sub-section 6.3)\n", + " 3. Compute normalization of the array (sub-section 6.4)\n", + " 4. Create the $\\mathbf{A}$ oracle operator for encoding the array (see sub-section 6.5)\n", + " 5. Create the correspondient Grover-like operator, $\\mathbf{G}(\\mathbf{A(f_{x_i})})$, from $\\mathbf{A}$ (see section 7)\n", + " 6. Execute the tested *amplitude estimation* algorithm properly confured using the operators $\\mathbf{G}(\\mathbf{A(f_{x_i})})$, and $\\mathbf{A}$ (see section 8)\n", + " 7. With the *amplitude estimation* algorithm result compute the integral (see section 8)\n", + " 8. Compute the desired metrics (see section 9).\n", + " 2. For each of the n qbits tested computed the mean, standard deviation, minimum and maximum of the 10 values of each metric.\n", + "2. For each one of the interval integration return the computed metric statistics for each number of qbits. \n", + "\n", + "\n", + "If a configuration parameter of the *amplitude estimation* algorithm want to be benchmarked for each posible values the complete above procedure should be executed.\n", + "\n", + "Generating a report using the different results of the benchmarking." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": "Python 3 (ipykernel)", + "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.11.7" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +}