Removed oracle_function as parameter for IQAE, adjusted QMCI, documen…

…tation, tests. Rework QMCI tutorial.

documentation/source/general/tutorial/QMCItutorial.rst

@@ -3,11 +3,11 @@
 Quantum Monte Carlo Integration with Iterative QAE
 ==================================================
 
-This tutorial will provide you with an introduction of Quantum Monte Carlo Integration within Qrisp.
+This tutorial will provide you with an introduction to Quantum Monte Carlo Integration within Qrisp.
 
-For this purpose we will first give you a theoretical overview of what this technique is about and where it is used. 
-Then we will dive into the practical implemention within Qrisp. This also includes the usage of :ref:`Iterative Quantum Amplitude estimation <IQAE>` . 
-To finish of this tutorial, we investigate the full implementation of simple example by integrating :math:`x^2` over a uniform distribution in the interval :math:`\lbrack 0,1 \rbrack` .
+For this purpose, we will first give you a theoretical overview of what this technique is about and where it is used. 
+Then we will dive into the practical implemention within Qrisp. This also includes the usage of :ref:`Iterative Quantum Amplitude Estimation <IQAE>`. 
+To finish of this tutorial, we investigate the full implementation of a simple example by integrating :math:`f(x)=x^2` over a uniform distribution in the interval :math:`\lbrack 0,1 \rbrack`.
 
 The relevant literature can be found in the following papers: `A general quantum algorithm for numerical integration <https://www.nature.com/articles/s41598-024-61010-9>`_ and `Option pricing using Quantum computers <https://arxiv.org/pdf/1905.02666>`_ for QMCI and `Accelerated Quantum Amplitude Estimation
 without QFT <https://arxiv.org/pdf/2407.16795>`_ for IQAE.
@@ -25,159 +25,90 @@ Mathemically speaking, we want to find an approximation for the following (gener
 
 As one does, we approximate an integral which cannot be solved analytically as a bunch of sums.
 
-Integrals like these may appear in many different places, from chemistry, over many-body physics to mathematical finance.
+Integrals like these may appear in many different places, from chemistry, through many-body physics, to mathematical finance.
 
 Implementation in Qrisp
 -----------------------
 
-The implementation in Qrisp requires the implementation of the function on a :ref:`QuantumFloat` . We will see how this can be done later in the example. 
+The implementation in Qrisp requires the implementation of the function on a :ref:`QuantumFloat`. We will see how this can be done later in the example. 
 
-There is multiple ways to implement Monte Carlo integration in a quantum fashion. Within Qrisp, we use the approach of :ref:`quantum counting <QCounting>` . The idea is to discretize not only the :math:`x`-axis but the :math:`y`-axis aswell. We use two ``QuantumFloats``  for this. 
-One ``QuantumFloat`` will hold the discretized values of the distribution, i.e. the relevant support on the :math:`x`-axis, while the other one will hold the discretized values of the :math:`y`-axis into which we will encode the the function values.
-We can then simply count, in quantum fashion, the number of points under the function curve, and divide it by the number of total points.
+There is multiple ways to implement Monte Carlo integration in a quantum fashion. Within Qrisp, we use an approach based on quantum counting. The idea is to discretize not only the :math:`x`-axis but the :math:`y`-axis as well. We use two ``QuantumFloats`` for this. 
+One ``QuantumFloat`` will hold the discretized values of the distribution, i.e., the relevant support on the :math:`x`-axis, while the other one will hold the discretized values of the :math:`y`-axis which we will encode the function values.
+We can then simply count, in quantum fashion, the number of points under the function curve, and divide it by the number of total points. This is facilitated by a ``QuantumBool`` which for each state $\ket{x}\ket{y}$ indicates whether the condition $y<f(x)$ is satisfied.
 
 Don't give up just yet, the mathematical description will bring you more clarity!
 
-We start with a ``state_function`` that encodes the ``QuantumFloats`` as follows
+We start with a ``state_function`` that, as a first step, encodes the ``QuantumFloats`` as follows:
 
 .. math::
 
-    \ket{0} \ket{0} \rightarrow \frac{1}{\sqrt{M \cdot N}} \sum^{N-1}_{x=0} \sum^{M-1}_{1=0} \ket{x} \ket{y}
+    \ket{0} \ket{0} \rightarrow \frac{1}{\sqrt{M \cdot N}} \sum^{N-1}_{x=0} \sum^{M-1}_{y=0} \ket{x} \ket{y}
 
-The ``oracle_function`` encodes the "number of points under the curve" condition as follows:
+As a second step, it applies the oracle that encodes the "points under the curve" condition as follows:
 
 .. math::
 
-    \ket{x} \ket{y} \rightarrow (-1)^{f(x) \leq y} \ket{x} \ket{y}
+    \ket{x} \ket{y} \ket{\text{False}} \rightarrow \mathbb{1}_{y \geq f(x)} \ket{x} \ket{y} \ket{\text{False}} + \mathbb{1}_{y < f(x)} \ket{x} \ket{y} \ket{\text{True}}
 
 We now arrive at the central step of this algorithm, which is :ref:`Quantum Amplitude Estimation <QAE>`. We use it to find
 
 .. math::
 
-    p(\{ (x,y) \mid f(x) \leq y \}) = \frac{1}{N} \sum^{N-1}_{x=0} \frac{1}{M} \sum^{M-1}_{y=0}  \mathbb{1}_{f(x) \leq y} \approx \frac{1}{N} \sum^{N-1}_{x=0} \frac{f(x)}{M}
+    p(\{ (x,y) \mid y < f(x) \}) = \frac{1}{N} \sum^{N-1}_{x=0} \frac{1}{M} \sum^{M-1}_{y=0}  \mathbb{1}_{y < f(x)} \approx \frac{1}{N} \sum^{N-1}_{x=0} \frac{f(x)}{M}
 
-The last expression is then the approximation for the integral in question. 
+The last expression is then (up to a scaling factor) the approximation for the integral in question. 
 
 
 Iterative Quantum Amplitude Estimation
 --------------------------------------
 
-In Qrisp we have the option of using an alternative amplitude estimation algorithm, namely `Accelerated Quantum Amplitude Estimation, see Algorithm 1 <https://arxiv.org/pdf/2407.16795>`_ , which iteratively applies :ref:`amplitude amplification <QAE>` to find an estimation for the amplitude of the good state.
+In Qrisp we have the option of using a resource efficient amplitude estimation algorithm, namely `Accelerated Quantum Amplitude Estimation, see Algorithm 1 <https://arxiv.org/pdf/2407.16795>`_ , which iteratively applies :ref:`amplitude amplification <QAE>` to find an estimation for the amplitude of the good state.
 The goal of the algorithm is as follows: 
 
-We start with a unitary operator :math:`\textbf{A}`, which represents a quantum circuit on :math:`w+1` qubits s.t.
+We start with a unitary operator :math:`\mathcal{A}`, which acts on the input QuantumVariables as
 
 .. math::
 
-    \textbf{A} \ket{0}_{w+1} = \sqrt{1-a} \ket{\psi_0}_{w} \ket{0} + \sqrt{a} \ket{\psi_1}_{w} \ket{1},
+    \textbf{A} \ket{0}\ket{\text{False}} = \sqrt{1-a} \ket{\Psi_0} \ket{\text{False}} + \sqrt{a} \ket{\Psi_1} \ket{\text{True}},
 
-where :math:`a \in \lbrack 0 , 1 \rbrack` is unknown. :math:`\ket{\psi_{0,1}}_{w}` are :math:`w`-qubit states. 
+where :math:`a \in [0,1]` is unknown.
 
 The algorithm in question allows for us to establish an estimate :math:`\hat{a}` of the unknown :math:`a`. 
 
-Mathematically speaking this means, given an error :math:`\epsilon` and a confidence level :math:`\alpha`, the accelerated Iterative Quantum Amplitude estimation finds an estimate :math:`\hat{a}`, s.t.
+Mathematically speaking this means, given an error :math:`\epsilon` and a confidence level :math:`\alpha`, the Accelerated Quantum Amplitude Estimation finds an estimate :math:`\hat{a}` such that
 
 .. math::
 
-    \text{P} \{ \mid \hat{a} - a \mid \leq \epsilon \} \geq 1 - \alpha 
+    \mathbb{P}\{|\hat{a} - a|\leq\epsilon\}\geq 1-\alpha 
 
-Below you can see the base function with implements this algorithm, including explanatory comments. It is a straight-forward translation from the theoretical ideas presented in the paper. For further explanations, have a look at the paper itself!
+A documentation explaining how to use the Qrisp implementation of this algorithm can found in the :ref:`IQAE <IQAE>` reference.
 
-The implementations of subroutines can found in the :ref:`accelerated IQAE <IQAE>` reference.
 
-::
-
-    def acc_IQAE(qargs,state_function, oracle_function, eps, alpha, kwargs_oracle = {}):
-        
-        # start by defining the relevant constants 
-        E = 1/2 * pow(np.sin(np.pi * 3/14), 2) -  1/2 * pow(np.sin(np.pi * 1/6), 2) 
-        F = 1/2 * np.arcsin(np.sqrt(2 * E))
-        C = 4/ (6*F + np.pi)
-
-        # the break condition defines when the algorithm converges with the desired accurarcy
-        break_cond =  2 * eps + 1
-        K_i = 1
-        m_i = 0
-        index_tot = 0
-        
-        # the main loop
-        while break_cond > 2 * eps : 
-            index_tot +=1
-            
-            # further constant defined
-            alp_i = C*alpha * eps * K_i 
-            N_i = int(np.ceil(1/(2 * pow(E, 2) ) * np.log(2/alp_i) ) )
-
-            # perform Quantum Amplitude amplification, and measure the number of |1> for the last qubit
-            qargs_dupl = [qarg.duplicate() for qarg in qargs]
-            A_i  = quantCirc( int((K_i -1 )/2) , N_i, qargs_dupl, state_function, 
-                            oracle_function, kwargs_oracle ) 
-            
-            for qarg in qargs_dupl:
-                qarg.delete()
-
-            
-            # compute new thetas
-            theta_b, theta_sh = compute_thetas(m_i,  K_i, A_i, E)
-            # compute new Li
-            L_new, m_new = compute_Li(m_i , K_i, theta_b, theta_sh)
-            
-            # assign new parameters
-            m_i = m_new
-            K_i = L_new * K_i
-            
-            # set new breaking condition
-            break_cond = abs( theta_b - theta_sh )
-        
-        # return the final approximation 
-        final_res = np.sin((theta_b+theta_sh)/2)**2
-        return final_res
-
-
-
-
-The QMCI class - full example
------------------------------
-
-Next up, we will go through a full example implementation to integrate :math:`x^2` over a uniform distribution in :math:`\lbrack 0,1 \rbrack`. This is the equivalent to the QMCI function. 
-
-First, we define the uniform distribution on a ``QuantumFloat``, which is just a uniform superposition of all qubits.
-
-::
-
-    def uniform(*args):
-        for arg in args:
-            h(arg)
-
-We also need a function that we want to integrate.
-
-::
-
-    def f(qf):
-        return qf*qf
-    
+QMCI example implementation
+---------------------------
 
-Next, we create the ``QuantumFloat``, on which we evaluate our function and a duplicate for the discretization of the :math:`y`-axis
+Next up, we will step-by-step go through a full example implementation of QMCI tailored to the example of integrating the function $f(x)=x^2$ over the uniform distribution in the interval $[0,1]$.
+A general implementation for integration of multidimensional functions is provided by the :ref:`QMCI <QMCI>` function.
 
+First, we define the function that we want to integrate, and a function for preparing the uniform distribution. 
+Additionally, we define a list of variables ``qargs`` repesenting the $x$-axis (``qargs[0]``) and $y$-axis (``qargs[1]``). 
+Thereby, the QuantumVariable representing the $y$-axis has to be chosen appropriately with respect to the values that the result of ``function(qargs[0])`` assumes.
 
 ::
 
-    
-    qf = QuantumFloat(2,-2)
+    from qrisp import *
 
-    dupl_args = [arg.duplicate() for arg in qargs]
-    dupl_res_qf = function(*dupl_args)
-    qargs.append(dupl_res_qf.duplicate())
+    def function(qf):
+        return qf*qf
 
-    for arg in dupl_args:
-        arg.delete()
-    dupl_res_qf.delete()
+    def distribution(qf):
+        h(qf)
 
+    qargs = [QuantumFloat(2,-2), QuantumFloat(4,-4)]
 
-We also have consider whether the ``QuantumFloat`` is not definded within a interval that differs from :math:`\lbrack 0, 1 \rbrack` . 
-In a way we calculate the volume of space over which the ``QuantumFloat`` is defined.
+Second, we determine the correct scaling factor by calculating the volume of the hypercube spanned by the intervals for the $x$-axis and $y$-axis.
 
-We also append a ``QuantumBool`` to our input ``qargs``, which will serve as the final qubit to be measured, i.e. the qubit in register :math:`w+1`.  
+We also append a ``QuantumBool`` to our input ``qargs``, which will indicate the "points under the curve". 
 
 ::
 
@@ -187,9 +118,7 @@ We also append a ``QuantumBool`` to our input ``qargs``, which will serve as the
     
     qargs.append(QuantumBool())
 
-Now we arrive at the heart of the algorithm, the definition of the ``oracle_function`` and the ``state_function``.
-
-Let us first look at the ``state_function``:
+Now, we arrive at the heart of the algorithm, the definition of the ``state_function``:
 
 ::
 
@@ -201,30 +130,23 @@ Let us first look at the ``state_function``:
 
         distribution(qf_x)
         h(qf_y)
+
         qbl = (qf_y < function(qf_x))
         cx(qbl,tar)
 
-It receives the ``@auto_uncompute`` :ref:`decorator <uncomputation>`. We apply the chosen distribution to ``qf_x``, which represents the :math:`x`-axis support. As explained earlier, we also discretize the :math:`y`-axis by appling an ``h``-gate to ``qf_y``.
-We then evaluate in superposition which states in ``qf_y`` are smaller than the chosen function acting on ``qf_x``, i.e. the function's support in the distribution.
+It receives the ``@auto_uncompute`` :ref:`decorator <uncomputation>` ensuring that all intermediate variables are properly uncomputed. 
+We apply the chosen distribution to ``qf_x``, which represents the :math:`x`-axes support. 
+As explained earlier, we also discretize the :math:`y`-axis by appling an ``h`` gate to ``qf_y``.
+We then evaluate in superposition which states in ``qf_y`` are smaller than the chosen function evaluated on ``qf_x``.
 
-We save the result of the comparison in a ``QuantumBool``, from which we can extract the measurement of the final qubit in register :math:`w+1` by applying a ``cx`` gate on the previously mentioned ``QuantumBool``
+We store the result of the comparison in the QuantumBool ``tar``, by applying a ``cx`` gate on the previously mentioned QuantumBool.
 
-This leads us to the ``oracle_function``
+With everything in place, we can now execute the :ref:`Iterative QAE algorithm <IQAE>`, with a chosen error tolerance ``eps`` and a confidence level ``alpha``.
+We also have to rescale with the previously calculated volume ``V0``.
 
 ::
-    
-    def oracle_function(*args):  
-        tar = args[2]
-        z(tar)
 
-It simply serves the function of tagging the :math:`\ket{1}`-state of the final qubit.
-
-With everything in place we can now execute the Iterative QAE algorithm, with a chosen error tolerance ``eps`` and a confidence level ``alpha``
-We also have to rescale with the previously calculated volume ``V0`` .
-
-::
-
-    a = acc_QAE(qargs, state_function, oracle_function, eps= 0.01, alpha= 0.01) 
+    a = IQAE(qargs, state_function, eps=0.01, alpha=0.01) 
     V = V0*a
 
 Aaaand that's it! The QMCI is complete! 
@@ -238,4 +160,3 @@ Let us now have a look at the result, and compare it to the expected result:
 
     >>> (0+0.25**2+0.5**2+0.75**2)/4
     0.21855991519015455
-

src/qrisp/alg_primitives/iterative_qae.py

@@ -16,40 +16,35 @@
 ********************************************************************************/
 """
 
+from qrisp import z
 from qrisp.alg_primitives.qae import amplitude_amplification 
 import numpy as np
 
 
-def IQAE(qargs, state_function, oracle_function, eps, alpha, kwargs_oracle = {}):
+def IQAE(qargs, state_function, eps, alpha, kwargs_oracle = {}):
     r"""
     Accelerated Quantum Amplitude Estimation (IQAE). This function performs :ref:`QAE <QAE>` with a fraction of the quantum resources of the well-known `QAE algorithm <https://arxiv.org/abs/quant-ph/0005055>`_.
     See `Accelerated Quantum Amplitude Estimation without QFT <https://arxiv.org/abs/2407.16795>`_.
 
-    The problem of quantum amplitude estimation is described as follows:
+    The problem of iterative quantum amplitude estimation is described as follows:
 
-    * Given a unitary operator :math:`\mathcal{A}`, let :math:`\ket{\Psi}=\mathcal{A}\ket{0}`.
-    * Write :math:`\ket{\Psi}=\ket{\Psi_1}+\ket{\Psi_0}` as a superposition of the orthogonal good and bad components of :math:`\ket{\Psi}`.
-    * Find an estimate for :math:`a=\langle\Psi_1|\Psi_1\rangle`, the probability that a measurement of $\ket{\Psi}$ yields a good state.
+    * Given a unitary operator :math:`\mathcal{A}`, let :math:`\ket{\Psi}=\mathcal{A}\ket{0}\ket{\text{False}}`.
+    * Write :math:`\ket{\Psi}=\sqrt{a}\ket{\Psi_1}\ket{\text{True}}+\sqrt{1-a}\ket{\Psi_0}\ket{\text{False}}` as a superposition of the orthogonal good and bad components of :math:`\ket{\Psi}`.
+    * Find an estimate for $a$, the probability that a measurement of $\ket{\Psi}$ yields a good state.
 
     Parameters
     ----------
-    qargs : QuantumVariable or list[QuantumVariable]
-        The (list of) QuantumVariables which represent the state,
-        the quantum amplitude estimation is performed on.
+    qargs : list[QuantumVariable]
+        The list of QuantumVariables which represent the state,
+        the quantum amplitude estimation is performed on. The last variable in the list must be of type :ref:`QuantumBool`.
     state_function : function
         A Python function preparing the state :math:`\ket{\Psi}`.
         This function will receive the variables in the list ``qargs`` as arguments in the
         course of this algorithm.
-    oracle_function : function
-        A Python function tagging the good state :math:`\ket{\Psi_1}`.
-        This function will receive the variables in the list ``qargs`` as arguments in the
-        course of this algorithm.
     eps : float
         Accuracy $\epsilon>0$ of the algorithm.
     alpha : float
         Confidence level $\alpha\in (0,1)$ of the algorithm.
-    kwargs_oracle : dict, optional
-        A dictionary containing keyword arguments for the oracle. The default is {}.
 
     Returns
     -------
@@ -95,27 +90,22 @@ def state_function(inp, tar):
                 with control(inp[k]):
                     ry(2**(k+1)/N,tar)
 
-    We define the ``oracle_function``:
-
-    ::
-
-        def oracle_function(inp, tar):
-            z(tar)
-
     Finally, we apply IQAE and obtain an estimate $a$ for the value of the integral $A=0.27268$.
 
     ::
 
-        eps = 0.01
-        alpha = 0.01
-
-        a = IQAE(input_list, state_function, oracle_function, eps=eps, alpha=alpha)
+        a = IQAE(input_list, state_function, eps=0.01, alpha=0.01)
 
     >>> a 
     0.26782038552705856
 
     """
 
+    # The oracle tagging the good states
+    def oracle_function(*args):  
+        tar = args[-1]
+        z(tar)
+
     E = 1/2 * pow(np.sin(np.pi * 3/14), 2) -  1/2 * pow(np.sin(np.pi * 1/6), 2) 
     F = 1/2 * np.arcsin(np.sqrt(2 * E))