From 6bc2bd9787b7b01d44c970cfd5a352155dd40821 Mon Sep 17 00:00:00 2001 From: gonfeco Date: Mon, 13 May 2024 20:36:00 +0200 Subject: [PATCH] Added updated modules from QQuantLib of the FinancialApoplications software package --- tnbs/BTC_02_AE/QQuantLib/AE/ae_class.py | 82 +- .../QQuantLib/AE/ae_classical_qpe.py | 11 +- .../QQuantLib/AE/ae_iterative_quantum_pe.py | 11 +- .../QQuantLib/AE/extended_real_quantum_ae.py | 898 ++++++++++++++++++ .../QQuantLib/AE/iterative_quantum_ae.py | 244 +++-- .../QQuantLib/AE/maximum_likelihood_ae.py | 23 +- .../AE/modified_iterative_quantum_ae.py | 459 +++++++++ .../QQuantLib/AE/modified_real_quantum_ae.py | 563 +++++++++++ tnbs/BTC_02_AE/QQuantLib/AE/montecarlo_ae.py | 17 +- .../BTC_02_AE/QQuantLib/AE/real_quantum_ae.py | 132 ++- .../QQuantLib/AE/shots_real_quantum_ae.py | 569 +++++++++++ tnbs/BTC_02_AE/QQuantLib/DL/data_loading.py | 16 +- .../QQuantLib/DL/encoding_protocols.py | 15 +- tnbs/BTC_02_AE/QQuantLib/PE/classical_qpe.py | 10 +- .../QQuantLib/PE/iterative_quantum_pe.py | 9 +- .../QQuantLib/finance/ae_price_estimation.py | 162 ---- .../QQuantLib/finance/classical_finance.py | 691 -------------- .../QQuantLib/finance/payoff_class.py | 101 -- .../QQuantLib/finance/probability_class.py | 144 --- .../QQuantLib/finance/quantum_integration.py | 49 +- .../01_AmplitudeEstimationKernel.ipynb | 10 +- .../notebooks/02_AmplitudeEstimationBTC.ipynb | 35 +- ...mplitudeEstimationBenchmarkExecution.ipynb | 6 +- .../04_AmpliutdeEstimationAlgorithms.ipynb | 20 +- .../QQuantLib/utils/data_extracting.py | 8 +- tnbs/BTC_02_AE/QQuantLib/utils/qlm_solver.py | 57 -- tnbs/BTC_02_AE/QQuantLib/utils/utils.py | 4 +- tnbs/BTC_02_AE/ae_sine_integral.py | 169 +++- 28 files changed, 3027 insertions(+), 1488 deletions(-) create mode 100644 tnbs/BTC_02_AE/QQuantLib/AE/extended_real_quantum_ae.py create mode 100644 tnbs/BTC_02_AE/QQuantLib/AE/modified_iterative_quantum_ae.py create mode 100644 tnbs/BTC_02_AE/QQuantLib/AE/modified_real_quantum_ae.py create mode 100644 tnbs/BTC_02_AE/QQuantLib/AE/shots_real_quantum_ae.py delete mode 100644 tnbs/BTC_02_AE/QQuantLib/finance/ae_price_estimation.py delete mode 100644 tnbs/BTC_02_AE/QQuantLib/finance/classical_finance.py delete mode 100644 tnbs/BTC_02_AE/QQuantLib/finance/payoff_class.py delete mode 100644 tnbs/BTC_02_AE/QQuantLib/finance/probability_class.py delete mode 100644 tnbs/BTC_02_AE/QQuantLib/utils/qlm_solver.py diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/ae_class.py b/tnbs/BTC_02_AE/QQuantLib/AE/ae_class.py index 6aa67ab..dee84f1 100644 --- a/tnbs/BTC_02_AE/QQuantLib/AE/ae_class.py +++ b/tnbs/BTC_02_AE/QQuantLib/AE/ae_class.py @@ -5,17 +5,16 @@ """ -import sys -import os -import re import pandas as pd - -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.AE.maximum_likelihood_ae import MLAE from QQuantLib.AE.ae_classical_qpe import CQPEAE from QQuantLib.AE.ae_iterative_quantum_pe import IQPEAE from QQuantLib.AE.iterative_quantum_ae import IQAE +from QQuantLib.AE.modified_iterative_quantum_ae import mIQAE from QQuantLib.AE.real_quantum_ae import RQAE +from QQuantLib.AE.shots_real_quantum_ae import sRQAE +from QQuantLib.AE.extended_real_quantum_ae import eRQAE +from QQuantLib.AE.modified_real_quantum_ae import mRQAE from QQuantLib.AE.montecarlo_ae import MCAE from QQuantLib.utils.utils import text_is_none @@ -40,7 +39,7 @@ class AE: dictionary that allows the configuration of the AE algorithm. \\ The different configration keys of the different AE algorithms \\ can be provided. - """ +""" def __init__(self, oracle=None, target=None, index=None, ae_type=None, **kwargs): """ @@ -57,20 +56,18 @@ def __init__(self, oracle=None, target=None, index=None, ae_type=None, **kwargs) self.index = index #Processing kwargs - self.kwargs = kwargs - self.linalg_qpu = self.kwargs.get("qpu", None) + self.solver_dict = kwargs + # self.kwargs = kwargs + self.linalg_qpu = self.solver_dict.get("qpu", None) - # Set the QPU to use - self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") self._ae_type = ae_type #attributes created self.solver_ae = None self.ae_pdf = None - self.solver_dict = None self.oracle_calls = None self.max_oracle_depth = None self.schedule_pdf = None @@ -93,7 +90,6 @@ def ae_type(self, stringvalue): self._ae_type = stringvalue self.solver_ae = None self.ae_pdf = None - self.solver_dict = None self.oracle_calls = None def create_ae_solver(self): @@ -102,17 +98,8 @@ def create_ae_solver(self): """ text_is_none(self.ae_type, "ae_type attribute", variable_type=str) #common ae settings - self.solver_dict = { - "mcz_qlm" : self.kwargs.get("mcz_qlm", True), - "qpu" : self.kwargs.get("qpu", None) - } - if self.ae_type == "MLAE": - for par in ["delta", "ns", "schedule"]: - val_par = self.kwargs.get(par) - if val_par is not None: - self.solver_dict.update({par : val_par}) self.solver_ae = MLAE( self.oracle, target=self.target, @@ -120,10 +107,6 @@ def create_ae_solver(self): **self.solver_dict ) elif self.ae_type == "CQPEAE": - for par in ["auxiliar_qbits_number", "shots"]: - val_par = self.kwargs.get(par) - if val_par is not None: - self.solver_dict.update({par : val_par}) self.solver_ae = CQPEAE( self.oracle, target=self.target, @@ -131,10 +114,6 @@ def create_ae_solver(self): **self.solver_dict ) elif self.ae_type == "IQPEAE": - for par in ["cbits_number", "shots"]: - val_par = self.kwargs.get(par) - if val_par is not None: - self.solver_dict.update({par : val_par}) self.solver_ae = IQPEAE( self.oracle, target=self.target, @@ -142,32 +121,48 @@ def create_ae_solver(self): **self.solver_dict ) elif self.ae_type == "IQAE": - for par in ["epsilon", "alpha", "shots"]: - val_par = self.kwargs.get(par) - if val_par is not None: - self.solver_dict.update({par : val_par}) self.solver_ae = IQAE( self.oracle, target=self.target, index=self.index, **self.solver_dict ) + elif self.ae_type == "mIQAE": + self.solver_ae = mIQAE( + self.oracle, + target=self.target, + index=self.index, + **self.solver_dict + ) elif self.ae_type == "RQAE": - for par in ["epsilon", "gamma", "q"]: - val_par = self.kwargs.get(par) - if val_par is not None: - self.solver_dict.update({par : val_par}) self.solver_ae = RQAE( self.oracle, target=self.target, index=self.index, **self.solver_dict ) + elif self.ae_type == "mRQAE": + self.solver_ae = mRQAE( + self.oracle, + target=self.target, + index=self.index, + **self.solver_dict + ) + elif self.ae_type == "eRQAE": + self.solver_ae = eRQAE( + self.oracle, + target=self.target, + index=self.index, + **self.solver_dict + ) + elif self.ae_type == "sRQAE": + self.solver_ae = sRQAE( + self.oracle, + target=self.target, + index=self.index, + **self.solver_dict + ) elif self.ae_type == "MCAE": - for par in ["shots"]: - val_par = self.kwargs.get(par) - if val_par is not None: - self.solver_dict.update({par : val_par}) self.solver_ae = MCAE( self.oracle, target=self.target, @@ -187,6 +182,7 @@ def run(self): self.oracle_calls = self.solver_ae.oracle_calls self.max_oracle_depth = self.solver_ae.max_oracle_depth self.schedule_pdf = self.solver_ae.schedule_pdf + # Recover amplitude estimation from ae_solver self.ae_pdf = pd.DataFrame( [self.solver_ae.ae, self.solver_ae.ae_l, self.solver_ae.ae_u], diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/ae_classical_qpe.py b/tnbs/BTC_02_AE/QQuantLib/AE/ae_classical_qpe.py index f710aa6..c173abb 100644 --- a/tnbs/BTC_02_AE/QQuantLib/AE/ae_classical_qpe.py +++ b/tnbs/BTC_02_AE/QQuantLib/AE/ae_classical_qpe.py @@ -16,14 +16,10 @@ """ -import sys -import os -import re import time +#from copy import deepcopy import numpy as np import qat.lang.AQASM as qlm - -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.PE.classical_qpe import CQPE from QQuantLib.AA.amplitude_amplification import grover from QQuantLib.utils.utils import check_list_type @@ -71,9 +67,9 @@ def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): # Set the QPU to use self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") self.auxiliar_qbits_number = kwargs.get("auxiliar_qbits_number", 8) self.shots = int(kwargs.get("shots", 100)) @@ -239,4 +235,5 @@ def run(self): self.max_oracle_depth = 2 ** (int(self.auxiliar_qbits_number)-1) + 1 self.quantum_times = self.cqpe.quantum_times self.quantum_time = sum(self.cqpe.quantum_times) + return self.ae diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/ae_iterative_quantum_pe.py b/tnbs/BTC_02_AE/QQuantLib/AE/ae_iterative_quantum_pe.py index d1f9ddf..c4ceb31 100644 --- a/tnbs/BTC_02_AE/QQuantLib/AE/ae_iterative_quantum_pe.py +++ b/tnbs/BTC_02_AE/QQuantLib/AE/ae_iterative_quantum_pe.py @@ -14,15 +14,10 @@ """ -import sys -import os -import re -import json import time +#from copy import deepcopy import numpy as np import qat.lang.AQASM as qlm - -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.PE.iterative_quantum_pe import IQPE from QQuantLib.AA.amplitude_amplification import grover from QQuantLib.utils.utils import check_list_type @@ -70,9 +65,9 @@ def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): # Set the QPU to use self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") self.cbits_number = kwargs.get("cbits_number", 8) self.shots = int(kwargs.get("shots", 100)) diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/extended_real_quantum_ae.py b/tnbs/BTC_02_AE/QQuantLib/AE/extended_real_quantum_ae.py new file mode 100644 index 0000000..783b551 --- /dev/null +++ b/tnbs/BTC_02_AE/QQuantLib/AE/extended_real_quantum_ae.py @@ -0,0 +1,898 @@ + +""" +This module contains necessary functions and classes to implement +extended Real Quantum Amplitude Estimation. This algorithm is an +extension of the RQAE paper: + + Manzano, A., Musso, D., Leitao, A. et al. + Real Quantum Amplitude Estimation + Preprint + +Author: Gonzalo Ferro Costas & Alberto Manzano Herrero + +""" +import time +#from copy import deepcopy +import numpy as np +import pandas as pd +import qat.lang.AQASM as qlm +from QQuantLib.AA.amplitude_amplification import grover +from QQuantLib.utils.data_extracting import get_results +from QQuantLib.utils.utils import measure_state_probability, bitfield_to_int, check_list_type, mask + + +def schedule_exponential_constant(epsilon=1.0e-2, gamma=0.05, ratio=2): + """ + Schedule for an exponential amplification schedule + and a constant gamma schedule + + Parameters + ---------- + epsilon: float + Desired error for the estimation + gamma : float + confidence level: it will be expected that the probability + of getting an error higher than epsilon will be lower than alpha + ratio : float + amplification ratio + + Returns + ---------- + + k_list : list + Schedule for the amplification at each step of the extended RQAE + (Grover operator applications) + gamma_list : list + Schedule for the confidence at each step of the extended RQAE + + """ + + # This is the maximum amplification of the algorithm + # it depends of the desired epsilon + bigk_max_next = 0.5 / np.arcsin(2 * epsilon) - 2.0 + k_max_next = np.ceil((bigk_max_next - 1) / 2) + + # The first step allways must be 0 because in the first step + # we do not amplify. We want the sign + k_ = 0 + big_k = 2 * k_ + 1 + # We want a exponential schedule for the amplification + k_next = ratio + bigk_next = 2 * k_next + 1 + k_list = [k_] + while bigk_next < bigk_max_next: + k_ = k_next + bigk = 2 * k_ + 1 + k_next = k_ * ratio + bigk_next = 2 * k_next + 1 + k_list.append(int(np.ceil(k_))) + # We add the maximum amplification to the schedule + k_list.append(k_max_next) + # For the gamma we want a constant schedule + gamma_i = gamma / len(k_list) + gamma_list = [gamma_i] * len(k_list) + return k_list, gamma_list + +def schedule_exponential_exponential(epsilon, gamma, ratio_epsilon, ratio_gamma): + """ + Schedule for an exponential amplification schedule + and a exponential gamma schedule + + Parameters + ---------- + epsilon: float + Desired error for the estimation + gamma : float + confidence level: it will be expected that the probability + of getting an error higher than epsilon will be lower than alpha + ratio_epsilon : float + amplification ratio (ratio for setting the k at each step). + Only positive ratios. + ratio_gamma : float + ratio for selecting the gamma at each step of the extended RQAE + ratios can be positive or negative + + Returns + ---------- + + k_list : list + Schedule for the amplification at each step of the extended RQAE + (Grover operator applications) + gamma_list : list + Schedule for the confidence at each step of the extended RQAE + + """ + # This is the maximum amplification of the algorithm + # it depends of the desired epsilon + bigk_max_next = 0.5 / np.arcsin(2 * epsilon) - 2.0 + k_max_next = np.ceil((bigk_max_next - 1) / 2) + # The first step allways must be 0 because in the first step + # we do not amplify. We want the sign + k_ = 0 + big_k = 2 * k_ + 1 + # We want a exponential schedule for the amplification + k_next = np.abs(ratio_epsilon) + bigk_next = 2 * k_next + 1 + k_list = [k_] + while bigk_next < bigk_max_next: + k_ = k_next + bigk = 2 * k_ + 1 + k_next = k_ * np.abs(ratio_epsilon) + bigk_next = 2 * k_next + 1 + k_list.append(k_) + #Schedule for gamma: we want exponential schedule + gamma_i = [np.abs(ratio_gamma) ** i for i in range(len(k_list))] + cte = np.sum(gamma_i) / gamma + gamma_list = [gm / cte for gm in gamma_i] + if ratio_gamma < 0: + gamma_list.reverse() + return k_list, gamma_list + +def schedule_linear_linear(epsilon, gamma, slope_epsilon, slope_gamma): + """ + Schedule for a lineal amplification schedule + and a lineal gamma schedule + + Parameters + ---------- + epsilon: float + Desired error for the estimation + gamma : float + confidence level: it will be expected that the probability + of getting an error higher than epsilon will be lower than alpha + slope_epsilon : float + amplification slope (slope for setting the k at each step). + Only positive slope. + slope_gamma : float + slope for selecting the gamma at each step of the extended RQAE + slope can be positive or negative + + Returns + ---------- + + k_list : list + Schedule for the amplification at each step of the extended RQAE + (Grover operator applications) + gamma_list : list + Schedule for the confidence at each step of the extended RQAE + """ + # This is the maximum amplification of the algorithm + # it depends of the desired epsilon + bigk_max_next = 0.5 / np.arcsin(2 * epsilon) - 2.0 + k_max_next = np.ceil((bigk_max_next - 1) / 2) + + # The first step allways must be 0 because in the first step + # we do not amplify. We want the sign + k_ = 0 + big_k = 2 * k_ + 1 + # We want a exponential schedule for the amplification + k_next = k_ + slope_epsilon + bigk_next = 2 * k_next + 1 + k_list = [k_] + + while bigk_next < bigk_max_next: + k_ = k_next + bigk = 2 * k_ + 1 + k_next = k_ + slope_epsilon + bigk_next = 2 * k_next + 1 + k_list.append(k_) + + k_list.append(k_max_next) + #Schedule for gamma: we want a linear schedule + gamma_i = [np.abs(slope_gamma) * i + 1 for i in range(len(k_list))] + cte = np.sum(gamma_i) / gamma + gamma_list = [gm / cte for gm in gamma_i] + if slope_gamma < 0: + gamma_list.reverse() + return k_list, gamma_list + +def schedule_linear_constant(epsilon, gamma, slope_epsilon): + """ + Schedule for a lineal amplification schedule + and a lineal gamma schedule + + Parameters + ---------- + epsilon: float + Desired error for the estimation + gamma : float + confidence level: it will be expected that the probability + of getting an error higher than epsilon will be lower than alpha + slope_epsilon : float + amplification slope (slope for setting the k at each step). + Only positive slope. + slope_gamma : float + slope for selecting the gamma at each step of the extended RQAE + slope can be positive or negative + + Returns + ---------- + + k_list : list + Schedule for the amplification at each step of the extended RQAE + (Grover operator applications) + gamma_list : list + Schedule for the confidence at each step of the extended RQAE + """ + # This is the maximum amplification of the algorithm + # it depends of the desired epsilon + bigk_max_next = 0.5 / np.arcsin(2 * epsilon) - 2.0 + k_max_next = np.ceil((bigk_max_next - 1) / 2) + + # The first step allways must be 0 because in the first step + # we do not amplify. We want the sign + k_ = 0 + big_k = 2 * k_ + 1 + # We want a exponential schedule for the amplification + k_next = k_ + slope_epsilon + bigk_next = 2 * k_next + 1 + k_list = [k_] + + while bigk_next < bigk_max_next: + k_ = k_next + bigk = 2 * k_ + 1 + k_next = k_ + slope_epsilon + bigk_next = 2 * k_next + 1 + k_list.append(k_) + + k_list.append(k_max_next) + # For the gamma we want a constant schedule + gamma_i = gamma / len(k_list) + gamma_list = [gamma_i] * len(k_list) + return k_list, gamma_list + +def select_schedule(erqae_schedule, epsilon, gamma): + """ + Scheduler selector. + Parameters + ---------- + erqae_schedule : dict + Dictionary with configuration for scheduler + epsilon: float + Desired error for the estimation + gamma : float + confidence level: it will be expected that the probability + of getting an error higher than epsilon will be lower than alpha + Returns + ---------- + + k_list : list + Schedule for the amplification at each step of the extended RQAE + (Grover operator applications) + gamma_list : list + Schedule for the confidence at each step of the extended RQAE + """ + if type(erqae_schedule) != dict: + raise TypeError("erqae_schedule MUST BE a dictionary") + + schedule_type = erqae_schedule["type"] + ratio_slope_k = erqae_schedule["ratio_slope_k"] + ratio_slope_gamma = erqae_schedule["ratio_slope_gamma"] + + if schedule_type == "exp_const": + k_list, gamma_list = schedule_exponential_constant( + epsilon, gamma, ratio_slope_k) + elif schedule_type == "exp_exp": + k_list, gamma_list = schedule_exponential_exponential( + epsilon, gamma, ratio_slope_k, ratio_slope_gamma) + elif schedule_type == "linear_linear": + k_list, gamma_list = schedule_linear_linear( + epsilon, gamma, ratio_slope_k, ratio_slope_gamma) + elif schedule_type == "linear_const": + k_list, gamma_list = schedule_linear_constant( + epsilon, gamma, ratio_slope_k) + else: + raise ValueError("Not valid schedule_type provided") + # k_list.append(k_list[-1]) + # gamma_list.append(gamma_list[-1]) + return k_list, gamma_list + + + +class eRQAE: + """ + Class for extended Real Quantum Amplitude Estimation (RQAE) + algorithm + + Parameters + ---------- + oracle: QLM gate + QLM gate with the Oracle for implementing the + Grover operator + target : list of ints + python list with the target for the amplitude estimation + index : list of ints + qubits which mark the register to do the amplitude + estimation + + kwars : dictionary + dictionary that allows the configuration of the IQAE algorithm: \\ + Implemented keys: + + qpu : QLM solver + solver for simulating the resulting circuits + epsilon : int + precision + gamma : float + accuracy + mcz_qlm : bool + for using or not QLM implementation of the multi controlled Z + gate + """ + + def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): + """ + + Method for initializing the class + + """ + ########################################### + # Setting attributes + self._oracle = oracle + self._target = check_list_type(target, int) + self._index = check_list_type(index, int) + + # Set the QPU to use + self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU + if self.linalg_qpu is None: + raise ValueError("Not QPU was provide. Please provide it!") + + self.epsilon = kwargs.get("epsilon", 0.01) + self.gamma = kwargs.get("gamma", 0.05) + self.mcz_qlm = kwargs.get("mcz_qlm", True) + self.save_circuits = kwargs.get("save_circuits", False) + + # Schedule + self.erqae_schedule = kwargs.get("erqae_schedule", None) + if self.erqae_schedule is None: + raise ValueError("erqae_schedule kwargs CAN NOT BE NONE") + + self.schedule_k, self.schedule_gamma = select_schedule( + self.erqae_schedule, self.epsilon, self.gamma) + # print(self.schedule_k) + # print(self.schedule_gamma) + + # Creating the grover operator + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + self.ae_l = None + self.ae_u = None + self.ae = None + self.circuit_statistics = None + self.time_pdf = None + self.run_time = None + self.schedule = {} + self.oracle_calls = None + self.max_oracle_depth = None + self.schedule_pdf = None + self.quantum_times = [] + self.quantum_time = None + self.circuit_dict = {} + self.info = None + + + ##################################################################### + @property + def oracle(self): + """ + creating oracle property + """ + return self._oracle + + @oracle.setter + def oracle(self, value): + """ + setter of the oracle property + """ + self._oracle = value + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def target(self): + """ + creating target property + """ + return self._target + + @target.setter + def target(self, value): + """ + setter of the target property + """ + self._target = check_list_type(value, int) + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def index(self): + """ + creating index property + """ + return self._index + + @index.setter + def index(self, value): + """ + setter of the index property + """ + self._index = check_list_type(value, int) + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def shifted_oracle(self): + """ + creating shifted_oracle property + """ + return self._shifted_oracle + + @shifted_oracle.setter + def shifted_oracle(self, shift): + """ + setter of the shifted_oracle property + + Parameters + ---------- + shift : float + shift for the oracle + """ + self._shifted_oracle = qlm.QRoutine() + wires = self._shifted_oracle.new_wires(self.oracle.arity + 1) + self._shifted_oracle.apply(qlm.H, wires[-1]) + self._shifted_oracle.apply( + qlm.RY(2 * np.arccos(shift)).ctrl(), wires[-1], wires[0] + ) + self._shifted_oracle.apply( + mask( + self.oracle.arity, + 2**self.oracle.arity - 1 - bitfield_to_int(self.target), + ).ctrl(), + wires[-1], + wires[: self.oracle.arity], + ) + self._shifted_oracle.apply(qlm.X, wires[-1]) + self._shifted_oracle.apply( + self._oracle.ctrl(), wires[-1], wires[: self._oracle.arity] + ) + self._shifted_oracle.apply(qlm.X, wires[-1]) + self._shifted_oracle.apply(qlm.H, wires[-1]) + + @staticmethod + def chebysev_bound(n_samples: int, gamma: float): + """ + Computes the length of the confidence interval for a given number + of samples n_samples and an accuracy gamma. + + Parameters + ---------- + n_samples : int + number of samples + gamma : float + accuracy + + Returns + ---------- + length of the confidence interval + """ + return np.sqrt(1 / (2 * n_samples) * np.log(2 / gamma)) + + def first_step(self, shift: float, shots: int, gamma: float): + """ + This function implements the first step of the RQAE paper. The result + is a first estimation of the desired amplitude. + + Parameters + ---------- + shift : float + shift for the first iteration + shots : int + number of measurements + gamma : float + accuracy + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + + self.shifted_oracle = 2 * shift + start = time.time() + results, circuit, _, _ = get_results( + self._shifted_oracle, self.linalg_qpu, shots=shots + ) + if self.save_circuits: + self.circuit_dict.update({"first_step": self._shifted_oracle}) + + end = time.time() + self.quantum_times.append(end-start) + start = time.time() + #probability_sum = results["Probability"].iloc[ + # bitfield_to_int([0] + list(self.target)) + #] + probability_sum = measure_state_probability( + results, [0] + list(self.target) + ) + + #probability_diff = results["Probability"].iloc[ + # bitfield_to_int([1] + list(self.target)) + #] + probability_diff = measure_state_probability( + results, [1] + list(self.target) + ) + epsilon_probability = eRQAE.chebysev_bound(shots, gamma) + + amplitude_max = np.minimum( + (probability_sum - probability_diff) / (4 * shift) + + epsilon_probability / (2 * np.abs(shift)), + 1.0, + ) + amplitude_min = np.maximum( + (probability_sum - probability_diff) / (4 * shift) + - epsilon_probability / (2 * np.abs(shift)), + -1.0, + ) + end = time.time() + first_step_time = end - start + + return [amplitude_min, amplitude_max], circuit + + def run_step(self, shift: float, shots: int, gamma: float, k: int): + """ + This function implements a step of the RQAE paper. The result + is a refined estimation of the desired amplitude. + + Parameters + ---------- + shift : float + shift for the first iteration + shots : int + number of measurements + gamma : float + accuracy + k : int + number of amplifications + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + #print(shift) + self.shifted_oracle = 2 * shift + + grover_oracle = grover( + self.shifted_oracle, + [0] + list(self.target), + np.arange(len(self.index) + 1), + mcz_qlm=self.mcz_qlm, + ) + + routine = qlm.QRoutine() + wires = routine.new_wires(self.shifted_oracle.arity) + routine.apply(self.shifted_oracle, wires) + for i in range(k): + routine.apply(grover_oracle, wires) + start = time.time() + results, circuit, _, _ = get_results(routine, self.linalg_qpu, shots=shots) + if self.save_circuits: + self.circuit_dict.update( + {"step_{}".format(k): routine} + ) + end = time.time() + self.quantum_times.append(end-start) + #probability_sum = results["Probability"].iloc[ + # bitfield_to_int([0] + list(self.target)) + #] + probability_sum = measure_state_probability( + results, [0] + list(self.target) + ) + + epsilon_probability = eRQAE.chebysev_bound(shots, gamma) + probability_max = min(probability_sum + epsilon_probability, 1) + probability_min = max(probability_sum - epsilon_probability, 0) + angle_max = np.arcsin(np.sqrt(probability_max)) / (2 * k + 1) + angle_min = np.arcsin(np.sqrt(probability_min)) / (2 * k + 1) + amplitude_max = np.sin(angle_max) - shift + amplitude_min = np.sin(angle_min) - shift + first_step_time = end - start + + return [amplitude_min, amplitude_max], circuit + + @staticmethod + def compute_info( + ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05, **kwargs + ): + """ + This function computes theoretical values of the IQAE algorithm. + + Parameters + ---------- + ratio: float + amplification ratio/policy + epsilon : float + precision + gamma : float + accuracy + + Return + ------ + info : dict + python dictionary with the computed information + + """ + epsilon = 0.5 * epsilon + theoretical_epsilon = None + k_max = None + big_t = None + gamma_i = None + n_grover = None + n_oracle = None + + info = { + "theoretical_epsilon": theoretical_epsilon, "k_max": k_max, + "big_t": big_t, "gamma_i": gamma_i, + "n_grover": n_grover, "n_oracle": n_oracle, + } + + return info + + @staticmethod + def display_information( + ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05, **kwargs + ): + """ + This function displays information of the properties of the + method for a given set of parameters + + Parameters + ---------- + ratio: float + amplification ratio/policy + epsilon : float + precision + gamma : float + accuracy + + """ + + info_dict = eRQAE.compute_info( + ratio=ratio, epsilon=epsilon, gamma=gamma) + + print("-------------------------------------------------------------") + print("BE AWARE: In extended RQAE the bounds depend on the selected schedule") + print("Here Not information is provided.") + print("Maximum number of amplifications: ", info_dict["k_max"]) + print("Maximum number of rounds: ", info_dict["big_t"]) + print("Maximum number of Grover operator calls: ", info_dict["n_grover"]) + print("Maximum number of Oracle operator calls: ", info_dict["n_oracle"]) + print("-------------------------------------------------------------") + + def erqae( self, epsilon: float = 0.01, gamma: float = 0.05, schedule_k: list = [], schedule_gamma: list = []): + """ + Implments the extended RQAE algorithm. The result is an estimation + of the desired amplitude with precision epsilon and accuracy gamma. + + Parameters + ---------- + epsilon : int + precision + gamma : float + accuracy + schedule_k : list + list with the amplification schedule + schedule_gamma : list + list with the confidence schedule + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + ###################################### + if len(schedule_k) == 0: + raise ValueError("The amplification schedule is empty!!") + if len(schedule_gamma) == 0: + raise ValueError("The gamma schedule is empty!!") + + epsilon = 0.5 * epsilon + # Always need to clean the circuit statistics property + self.circuit_statistics = {} + # time_list = [] + ############### First Step ####################### + shift = 0.5 # shift [-0.5, 0.5] + # For computing the desired epsilon to achieve we need the + # next step scheduled amplification + k_next = schedule_k[1] + bigk_next = 2 * k_next + 1 + # Probability epsilon for first step based on the next step + + # scheduled amplification + epsilon_first_p = np.abs(shift) * np.sin(0.5 * np.pi / (bigk_next + 2)) + + # Maximum probability epsilon for first step + epsilon_first_p_min = 0.5 * np.sin(0.5 * np.arcsin(2 * epsilon)) ** 2 + #print("epsilon: ", epsilon_first_p, "epsilon min: ", epsilon_first_p_min) + + #print(epsilon_first_p, epsilon_first_p_min, k_next) + # Real probabiliy epsilon to use for first step + epsilon_first_p = max(epsilon_first_p, epsilon_first_p_min) + # Gamma used for first step + gamma_first = schedule_gamma[0] + # This is shots for each iteration: Ni in the paper + n_first = int( + np.ceil(np.log(2 / gamma_first) / (2 * epsilon_first_p**2)) + ) + # Quantum routine for first step + [amplitude_min, amplitude_max], _ = self.first_step( + shift=shift, shots=n_first, gamma=gamma_first + ) + # print( + # "first step: epsilon_first_p: ", epsilon_first_p, + # "gamma_first: ", gamma_first, "n_first: ", n_first + # ) + self.schedule.update({0 : n_first}) + + # Real amplitude epsilon + epsilon_amplitude = (amplitude_max - amplitude_min) / 2 + + ############### Consecutive Steps ####################### + #print("i: ", 0, "epsilon_amplitude: ", epsilon_amplitude, "epsilon: ", epsilon) + i = 1 + while (epsilon_amplitude > epsilon) and (i < len(schedule_k) - 1): + + # This is the amplification for the current step + k_exp = int(np.floor(np.pi / (4 * np.arcsin(2 * epsilon_amplitude)) - 0.5)) + bigk_exp = 2 * k_exp + 1 + # Computation of the number of shots for the current iterative step + + # We need the schedule amplification of the current step and + # the schedule amplification of the following step + k_next = schedule_k[i+1] + bigk_next = 2 * k_next + 1 + k_current = schedule_k[i] + bigk_current = 2 * k_current + 1 + # We compute the desired epsilon for the current step + epsilon_step_p = 0.5 * np.sin( + 0.25 * np.pi * bigk_current / (bigk_next + 2)) ** 2 + # We compute the maximum epsilon achievable in the step + epsilon_step_p_min = 0.5 * np.sin( + 0.5 * bigk_exp * np.arcsin(2 * epsilon)) ** 2 + #print("epsilon: ", epsilon_step_p, "epsilon min: ", epsilon_step_p_min) + # Real probabiliy epsilon that should be achieved in the step + epsilon_step_p = max(epsilon_step_p, epsilon_step_p_min) + # gamma used for first step + gamma_step = schedule_gamma[i] + # This is the mandatory number of shots of the step for + # achieving the gamma_step and the epsilon_step_pT + n_step = int( + np.ceil(1 / (2 * epsilon_step_p**2) * np.log(2 / gamma_step)) + ) + if bigk_exp < bigk_current: + raise ValueError("bigk_exp: {} is lower than: bigk_current: {}".format( + bigk_exp, bigk_current) + ) + + # Quantum routine for current step + shift = -amplitude_min + if shift > 0: + shift = min(shift, 0.5) + if shift < 0: + shift = max(shift, -0.5) + # print("step: epsilon_step_p: ", epsilon_step_p, "gamma_step: ", + # gamma_step, "n_step: ", n_step, "k_exp: ", k_exp + # ) + [amplitude_min, amplitude_max], _ = self.run_step( + shift=shift, shots=n_step, gamma=gamma_step, k=k_exp + ) + # Added the shots for the k + if k_exp not in self.schedule: + self.schedule.update({k_exp:n_step}) + else: + # If k exists sum the shots with the before number of shots + self.schedule.update({k_exp:self.schedule[k_exp] + n_step}) + # time_list.append(time_pdf) + epsilon_amplitude = (amplitude_max - amplitude_min) / 2 + #print("i: ", i, "epsilon_amplitude: ", epsilon_amplitude, "epsilon: ", epsilon) + i = i + 1 + + if epsilon_amplitude > epsilon: + # This is the amplification for the current step + k_exp = int(np.floor(np.pi / (4 * np.arcsin(2 * epsilon_amplitude)) - 0.5)) + k_current = schedule_k[i] + bigk_current = 2 * k_current + 1 + bigk_exp = 2 * k_exp + 1 + # Computation of the number of shots for the current iterative step + + # We compute the maximum epsilon achievable in the step + epsilon_step_p = 0.5 * np.sin( + 0.5 * bigk_exp * np.arcsin(2 * epsilon)) ** 2 + #print("epsilon: ", epsilon_step_p, "epsilon min: ", epsilon_step_p_min) + # gamma used for first step + gamma_step = schedule_gamma[i] + # This is the mandatory number of shots of the step for + # achieving the gamma_step and the epsilon_step_pT + n_step = int( + np.ceil(1 / (2 * epsilon_step_p**2) * np.log(2 / gamma_step)) + ) + if bigk_exp < bigk_current: + raise ValueError("bigk_exp: {} is lower than: bigk_current: {}".format( + bigk_exp, bigk_current) + ) + + # Quantum routine for current step + shift = -amplitude_min + if shift > 0: + shift = min(shift, 0.5) + if shift < 0: + shift = max(shift, -0.5) + #print("step: epsilon_step_p: ", epsilon_step_p, "gamma_step: ", + # gamma_step, "n_step: ", n_step, "k_exp: ", k_exp + #) + [amplitude_min, amplitude_max], _ = self.run_step( + shift=shift, shots=n_step, gamma=gamma_step, k=k_exp + ) + # Added the shots for the k + if k_exp not in self.schedule: + self.schedule.update({k_exp:n_step}) + else: + # If k exists sum the shots with the before number of shots + self.schedule.update({k_exp:self.schedule[k_exp] + n_step}) + # time_list.append(time_pdf) + epsilon_amplitude = (amplitude_max - amplitude_min) / 2 + #print("i: ", i, "epsilon_amplitude: ", epsilon_amplitude, "epsilon: ", epsilon) + i = i + 1 + + return [2 * amplitude_min, 2 * amplitude_max] + + def run(self): + r""" + run method for the class. + + Returns + ---------- + + self.ae : + amplitude estimation parameter + + """ + start = time.time() + if len(self.schedule_k) == 0: + raise ValueError("The amplification schedule is empty!!") + if len(self.schedule_gamma) == 0: + raise ValueError("The gamma schedule is empty!!") + [self.ae_l, self.ae_u] = self.erqae( + epsilon=self.epsilon, gamma=self.gamma, + schedule_k=self.schedule_k, schedule_gamma=self.schedule_gamma + ) + self.ae = (self.ae_u + self.ae_l) / 2.0 + end = time.time() + self.run_time = end - start + self.schedule_pdf = pd.DataFrame.from_dict( + self.schedule, + columns=['shots'], + orient='index' + ) + self.schedule_pdf.reset_index(inplace=True) + self.schedule_pdf.rename(columns={'index': 'm_k'}, inplace=True) + self.oracle_calls = np.sum( + self.schedule_pdf['shots'] * (2 * self.schedule_pdf['m_k'] + 1)) + self.max_oracle_depth = np.max(2 * self.schedule_pdf['m_k']+ 1) + self.quantum_time = sum(self.quantum_times) + return self.ae diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/iterative_quantum_ae.py b/tnbs/BTC_02_AE/QQuantLib/AE/iterative_quantum_ae.py index fc488ea..eeb4b7f 100644 --- a/tnbs/BTC_02_AE/QQuantLib/AE/iterative_quantum_ae.py +++ b/tnbs/BTC_02_AE/QQuantLib/AE/iterative_quantum_ae.py @@ -11,20 +11,16 @@ """ -import sys -import os -import re import time +#from copy import deepcopy import numpy as np import pandas as pd import qat.lang.AQASM as qlm - - -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.AA.amplitude_amplification import grover from QQuantLib.utils.data_extracting import get_results -from QQuantLib.utils.utils import check_list_type, \ - measure_state_probability +from QQuantLib.utils.utils import check_list_type, measure_state_probability +import logging +logger = logging.getLogger('__name__') class IQAE: @@ -72,9 +68,9 @@ def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): # Set the QPU to use self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") self.epsilon = kwargs.get("epsilon", 0.01) self.alpha = kwargs.get("alpha", 0.05) self.shots = int(kwargs.get("shots", 100)) @@ -256,6 +252,50 @@ def invert_sector(a_min: float, a_max: float, flag: bool = True): return [theta_min, theta_max] + @staticmethod + def compute_info( + epsilon: float = 0.01, shots: int = 100, alpha: float = 0.05 + ): + """ + This function computes theoretical values of the IQAE algorithm. + + Parameters + ---------- + epsilon : float + precision + alpha : float + accuracy + shots : int + number of measurements on each iteration + + Return + ------ + info : dict + python dictionary with the computed information + + """ + # Maximum number of rounds: T in paper + big_t = np.ceil(np.log2(np.pi / (8 * epsilon))) + # Total number of Grover operator calls + n_grover = int( + 50 * np.log(2 / alpha * np.log2(np.pi / (4 * epsilon))) / epsilon + ) + # Maximum number of shots at each round + n_max = int(32 * np.log(2 / alpha * np.log2(np.pi / (4 * epsilon))) \ + / (1 - 2 * np.sin(np.pi / 14)) ** 2) + # This is the number of calls to the oracle operator (A) + n_oracle = 2 * n_grover + n_max * (big_t + 1) + # This is L in the papper + big_l = (np.arcsin(2 / shots * np.log(2 * big_t / epsilon))) ** 0.25 + k_max = big_l / epsilon / 2 + + info = { + "big_t": big_t, "n_grover": n_grover, "n_max": n_max, + "n_oracle": n_oracle, "big_l": big_l, "k_max": k_max + } + + return info + @staticmethod def display_information( epsilon: float = 0.01, shots: int = 100, alpha: float = 0.05 @@ -281,23 +321,14 @@ def display_information( print("N: ", shots) print("-------------------------------------------------------------") - # This is T, number of rounds, in the papper - big_t = np.ceil(np.log2(np.pi / (8 * epsilon))) - n_max = ( - 32 - / (1 - 2 * np.sin(np.pi / 14)) ** 2 - * np.log(2 / alpha * np.log2(np.pi / (4 * epsilon))) - ) - n_oracle = 50 / epsilon * np.log(2 / alpha * np.log2(np.pi / (4 * epsilon))) - # This is L in the papper - big_l = (np.arcsin(2 / shots * np.log(2 * big_t / epsilon))) ** 0.25 - k_max = big_l / epsilon / 2 + info_dict = IQAE.compute_info(epsilon, shots, alpha) print("-------------------------------------------------------------") - print("Maximum number of rounds: ", int(big_t)) - print("Maximum number of shots per round needed: ", int(n_max)) - print("Maximum number of amplifications: ", int(k_max)) - print("Maximum number of calls to the oracle: ", int(n_oracle)) + print("Maximum number of rounds: ", info_dict["big_t"]) + print("Maximum number of shots per round needed: ", info_dict["n_max"]) + print("Maximum number of amplifications: ", info_dict["k_max"]) + print("Maximum number of Grover operator calls: ", info_dict["n_grover"]) + print("Maximum number of Oracle operator calls: ", info_dict["n_oracle"]) print("-------------------------------------------------------------") @staticmethod @@ -355,78 +386,123 @@ def iqae(self, epsilon: float = 0.01, shots: int = 100, alpha: float = 0.05): [theta_l, theta_u] = [0.0, np.pi / 2] # This is T the number of rounds in the paper big_t = int(np.ceil(np.log2(np.pi / (8 * epsilon))) + 1) - # This is L in the paper - big_l = (np.arcsin(2 / shots * np.log(2 * big_t / epsilon))) ** 0.25 ##################################################### h_k = 0 n_effective = 0 # time_list = [] j = 0 # pure counter - - while theta_u - theta_l > 2 * epsilon: + # For avoiding infinite loop + failure_test = True + while (theta_u - theta_l > 2 * epsilon) and (failure_test == True): #start = time.time() i = i + 1 k_old = k [k, flag] = self.find_next_k(k_old, theta_l, theta_u, flag) big_k = 4 * k + 2 - #end = time.time() - #finding_time = end - start - - ##################################################### - routine = self.quantum_step(k) - start = time.time() - results, circuit, _, _ = get_results( - routine, linalg_qpu=self.linalg_qpu, shots=shots, qubits=self.index - ) - end = time.time() - self.quantum_times.append(end-start) - # time_pdf["m_k"] = k - a_ = measure_state_probability(results, self.target) - #a_ = results["Probability"].iloc[bitfield_to_int(self.target)] - ##################################################### - # Aggregate results from different iterations - if j == 0: - # In the first step we need to store the circuit statistics - step_circuit_stats = circuit.statistics() - step_circuit_stats.update({"n_shots": shots}) + + # The IQAE algorithm has a failure that can lead to infinite + # loops. In order to avoid it we are going to establish + # the following condition for executing it. The number of + # revolution of the 2 angles when apply k MUST BE the same + + # Number of revolutions that the 2 angles will execute + rounds_theta_l = np.floor(big_k * theta_l / (2 * np.pi)) + rounds_theta_u = np.floor(big_k * theta_u / (2 * np.pi)) + + if int(rounds_theta_l) == int(rounds_theta_u): + # If they are the same then we can execute the algorithm + ##################################################### + # Quantum Routine and Measurement + routine = self.quantum_step(k) + start = time.time() + results, circuit, _, _ = get_results( + routine, linalg_qpu=self.linalg_qpu, shots=shots, qubits=self.index + ) + end = time.time() + self.quantum_times.append(end-start) + # time_pdf["m_k"] = k + a_ = measure_state_probability(results, self.target) + ##################################################### + if j == 0: + # In the first step we need to store the circuit statistics + step_circuit_stats = circuit.statistics() + step_circuit_stats.update({"n_shots": shots}) + self.circuit_statistics.update({k: step_circuit_stats}) + self.schedule.update({k:shots}) + + # Aggregate results from different iterations + if k == k_old: + h_k = h_k + int(a_ * shots) + n_effective = n_effective + shots + a_ = h_k / n_effective + i = i - 1 + # Only update shots for the k application + step_circuit_stats = self.circuit_statistics[k] + step_circuit_stats.update({"n_shots": n_effective}) + self.schedule.update({k:n_effective}) + + else: + h_k = int(a_ * shots) + n_effective = shots + # Store the circuit statistics for new k + step_circuit_stats = circuit.statistics() + step_circuit_stats.update({"n_shots": shots}) + self.schedule.update({k:shots}) self.circuit_statistics.update({k: step_circuit_stats}) - self.schedule.update({k:shots}) - - if k == k_old: - h_k = h_k + int(a_ * shots) - n_effective = n_effective + shots - a_ = h_k / n_effective - i = i - 1 - # Only update shots for the k application - step_circuit_stats = self.circuit_statistics[k] - step_circuit_stats.update({"n_shots": n_effective}) - self.schedule.update({k:n_effective}) + # Compute the rest + epsilon_a = IQAE.chebysev_bound(n_effective, alpha / big_t) + a_max = np.minimum(a_ + epsilon_a, 1.0) + a_min = np.maximum(a_ - epsilon_a, 0.0) + + [theta_min, theta_max] = self.invert_sector(a_min, a_max, flag) + + #msg_txt = """Step j= {}. theta_l_before= {}. theta_u_before= + # k= {}, n_effective={}""".format( + # j, theta_l, theta_u, k, n_effective + # ) + #logger.warning(msg_txt) + + theta_l_ = ( + 2 * np.pi * np.floor(big_k * theta_l / (2 * np.pi)) + theta_min + ) / big_k + + theta_u_ = ( + 2 * np.pi * np.floor(big_k * theta_u / (2 * np.pi)) + theta_max + ) / big_k + + #msg_txt = """Step j= {}. theta_l_new= {}. + # theta_u_new= {}""".format( + # j, theta_l_, theta_u_ + # ) + #logger.warning(msg_txt) + + # If bounded limits are worse than step before limits use these ones + #theta_l = np.maximum(theta_l, theta_l_) + #theta_u = np.minimum(theta_u, theta_u_) + + #if (theta_l_ < theta_l) or (theta_u_ > theta_u): + # logger.warning("PROBLEM") + # msg_txt = """|_K * theta_l_|= {}. + # |_K * theta_u_|= {}""".format( + # np.floor(big_k * theta_l_ / (2 * np.pi)), + # np.floor(big_k * theta_u_ / (2 * np.pi)) + # ) + # logger.warning(msg_txt) + theta_l = theta_l_ + theta_u = theta_u_ + j = j + 1 else: - h_k = int(a_ * shots) - n_effective = shots - # Store the circuit statistics for new k - step_circuit_stats = circuit.statistics() - step_circuit_stats.update({"n_shots": shots}) - self.schedule.update({k:shots}) - self.circuit_statistics.update({k: step_circuit_stats}) - - # Compute the rest - epsilon_a = IQAE.chebysev_bound(n_effective, alpha / big_t) - a_max = np.minimum(a_ + epsilon_a, 1.0) - a_min = np.maximum(a_ - epsilon_a, 0.0) - [theta_min, theta_max] = self.invert_sector(a_min, a_max, flag) - - theta_l_ = ( - 2 * np.pi * np.floor(big_k * theta_l / (2 * np.pi)) + theta_min - ) / big_k - theta_u_ = ( - 2 * np.pi * np.floor(big_k * theta_u / (2 * np.pi)) + theta_max - ) / big_k - # If bounded limits are worse than step before limits use these ones - theta_l = np.maximum(theta_l, theta_l_) - theta_u = np.minimum(theta_u, theta_u_) - j = j + 1 + # In this case algorithm will begin an infintie loop + # logger.warning("PROBLEM-2") + # msg_txt = """k= {} K= {}. Revolutions for theta_l: {}. + # Revolutions for theta_u: {}""".format( + # k, big_k, rounds_theta_l, rounds_theta_u) + # logger.warning(msg_txt) + # msg_txt = "Epsilon: {}".format(theta_u - theta_l) + # logger.warning(msg_txt) + # logger.warning("\n") + failure_test = False [a_l, a_u] = [np.sin(theta_l) ** 2, np.sin(theta_u) ** 2] return [a_l, a_u] diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/maximum_likelihood_ae.py b/tnbs/BTC_02_AE/QQuantLib/AE/maximum_likelihood_ae.py index e53f348..2d1ae37 100644 --- a/tnbs/BTC_02_AE/QQuantLib/AE/maximum_likelihood_ae.py +++ b/tnbs/BTC_02_AE/QQuantLib/AE/maximum_likelihood_ae.py @@ -11,26 +11,17 @@ """ -import sys -import os import time -import re +#from copy import deepcopy from functools import partial import numpy as np import pandas as pd import scipy.optimize as so import qat.lang.AQASM as qlm - - - - -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.AA.amplitude_amplification import grover from QQuantLib.utils.data_extracting import get_results -from QQuantLib.utils.utils import measure_state_probability, \ - check_list_type, load_qn_gate -from QQuantLib.AE.mlae_utils import likelihood, log_likelihood, \ - cost_function +from QQuantLib.utils.utils import measure_state_probability, check_list_type, load_qn_gate +from QQuantLib.AE.mlae_utils import likelihood, log_likelihood, cost_function class MLAE: """ @@ -80,9 +71,9 @@ def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): # Set the QPU to use self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") ##delta for avoid problems in 0 and pi/2 theta limits self.delta = kwargs.get("delta", 1.0e-6) # ns for the brute force optimizer @@ -485,6 +476,10 @@ def run(self) -> float: ) self.theta = self.theta[0] self.ae = np.sin(self.theta) ** 2 + # MLAE does not provide upper and lower bounds. We fix to the + # estimated value + self.ae_l = self.ae + self.ae_u = self.ae result = self.ae end = time.time() self.run_time = end - start diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/modified_iterative_quantum_ae.py b/tnbs/BTC_02_AE/QQuantLib/AE/modified_iterative_quantum_ae.py new file mode 100644 index 0000000..36fbfd1 --- /dev/null +++ b/tnbs/BTC_02_AE/QQuantLib/AE/modified_iterative_quantum_ae.py @@ -0,0 +1,459 @@ +""" +This module contains necessary functions and classes to implement +Iterative Quantum Amplitude Estimation based on the paper: + + Grinko, D., Gacon, J., Zoufal, C. et al. + Iterative Quantum Amplitude Estimation + npj Quantum Inf 7, 52 (2021). + https://doi.org/10.1038/s41534-021-00379-1 + +Author: Gonzalo Ferro Costas & Alberto Manzano Herrero + +""" + +import time +#from copy import deepcopy +import numpy as np +import pandas as pd +import qat.lang.AQASM as qlm +from QQuantLib.AA.amplitude_amplification import grover +from QQuantLib.utils.data_extracting import get_results +from QQuantLib.utils.utils import check_list_type, measure_state_probability + + +class mIQAE: + """ + Class for Iterative Quantum Amplitude Estimation (IQAE) + algorithm + + Parameters + ---------- + oracle: QLM gate + QLM gate with the Oracle for implementing the + Grover operator + target : list of ints + python list with the target for the amplitude estimation + index : list of ints + qubits which mark the register to do the amplitude + estimation + + kwars : dictionary + dictionary that allows the configuration of the IQAE algorithm: \\ + Implemented keys: + + qpu : QLM solver + solver for simulating the resulting circuits + epsilon : float + precision + alpha : float + accuracy + shots : int + number of measurements on each iteration + mcz_qlm : bool + for using or not QLM implementation of the multi controlled Z + gate + """ + + def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): + """ + + Method for initializing the class + """ + # Setting attributes + self._oracle = oracle + self._target = check_list_type(target, int) + self._index = check_list_type(index, int) + + # Set the QPU to use + self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU + if self.linalg_qpu is None: + raise ValueError("Not QPU was provide. Please provide it!") + self.epsilon = kwargs.get("epsilon", 0.01) + self.alpha = kwargs.get("alpha", 0.05) + self.shots = int(kwargs.get("shots", 100)) + self.mcz_qlm = kwargs.get("mcz_qlm", True) + + # Creating the grover operator + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + self.ae_l = None + self.ae_u = None + self.theta_l = None + self.theta_u = None + self.theta = None + self.ae = None + self.circuit_statistics = None + self.time_pdf = None + self.run_time = None + self.schedule = {} + self.oracle_calls = None + self.max_oracle_depth = None + self.schedule_pdf = None + self.quantum_times = [] + self.quantum_time = None + + ##################################################################### + @property + def oracle(self): + """ + creating oracle property + """ + return self._oracle + + @oracle.setter + def oracle(self, value): + """ + setter of the oracle property + """ + self._oracle = value + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def target(self): + """ + creating target property + """ + return self._target + + @target.setter + def target(self, value): + """ + setter of the target property + """ + self._target = check_list_type(value, int) + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def index(self): + """ + creating index property + """ + return self._index + + @index.setter + def index(self, value): + """ + setter of the index property + """ + self._index = check_list_type(value, int) + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + ##################################################################### + + @staticmethod + def find_next_k( k: int, theta_lower: float, theta_upper: float): + """ + This is an implementation of Algorithm 2 from the mIQAE paper. + This function computes the next suitable k. + + Parameters + ---------- + k : int + number of times to apply the grover operator to the quantum circuit + theta_lower : float + lower bound for the estimation of the angle + theta_upper : float + upper bound for the estimation of the angle + + Returns + ---------- + k : int + number of times to apply the grover operator to the quantum circuit + """ + # This is K_i in the paper + bigk_i = 2 * k + 1 + # This is K in Algorithm 2 + big_k = np.floor(0.5 * np.pi / (theta_upper - theta_lower)) + + if big_k % 2 == 0: + #if K is even + big_k = big_k -1 + + while big_k >= 3.0 * bigk_i: + amplified_lower = np.floor(big_k * theta_lower / (0.5 * np.pi)) + amplified_upper = np.ceil(big_k * theta_upper / (0.5 * np.pi)) + if amplified_lower == (amplified_upper -1): + k_i = (big_k - 1) / 2.0 + return int(k_i) + big_k = big_k -2 + return int(k) + + @staticmethod + def compute_info( + epsilon: float = 0.01, shots: int = 100, alpha: float = 0.05 + ): + """ + This function computes theoretical values of the IQAE algorithm. + + Parameters + ---------- + epsilon : float + precision + alpha : float + accuracy + shots : int + number of measurements on each iteration + + Return + ------ + info : dict + python dictionary with the computed information + + """ + + + # Upper bound for amplification: Kmax in paper + bigk_max = int(np.pi / (4.0 * epsilon)) + # Maximum number of rounds: log3(pi/ (4 * epsilon)) in paper + big_t = int(np.ceil(np.log2(np.pi / (4 * epsilon)) / np.log2(3))) + # constant C in the paper + big_c = 1.0 / ((np.sin(np.pi / 21.0) * np.sin(8.0 * np.pi /21.0)) ** 2) + # Total number of Grover operator calls: 3/2 * C * K_max * ln(sqrt(27)/alpa) + n_grover = int(1.5 * bigk_max * big_c * np.log(np.sqrt(27) / alpha)) + # Total number of oracle operator calls: 2 * n_grover + N0_max + # We need to compute maximum number of shots for k=0 + big_k_0 = 1 + alpha_0 = 2 * alpha * big_k_0 / (3 * bigk_max) + big_n_0 = 2.0 * big_c * np.log(2.0 / alpha_0) + # c1 = np.log(3 * bigk_max / alpha) + # c2 = 0.5 * (big_t + 1) * big_t * np.log(3) + # c3 = big_t * np.log(alpha) + n_oracle = 2.0 * n_grover + big_n_0 + + + info = { + "bigk_max": bigk_max, "big_t": big_t, "n_grover": n_grover, + "n_oracle": n_oracle, "big_n_0": big_n_0 + } + + return info + + @staticmethod + def display_information( + epsilon: float = 0.01, shots: int = 100, alpha: float = 0.05 + ): + """ + This function displays information of the properties of the + method for a given set of parameters + + Parameters + ---------- + epsilon : float + precision + alpha : float + accuracy + shots : int + number of measurements on each iteration + + """ + + print("-------------------------------------------------------------") + print("epsilon: ", epsilon) + print("alpha: ", alpha) + print("N: ", shots) + print("-------------------------------------------------------------") + + info_dict = mIQAE.compute_info(epsilon, shots, alpha) + + print("-------------------------------------------------------------") + print("Maximum amplification (Kmax)", info_dict["bigk_max"]) + print("Maximum number of rounds: ", info_dict["big_t"]) + print("Maximum number of Grover operator calls: ", info_dict["n_grover"]) + print("Maximum number of Oracle operator calls: ", info_dict["n_oracle"]) + print("-------------------------------------------------------------") + + @staticmethod + def chebysev_bound(n_samples: int, gamma: float): + r""" + Computes the length of the confidence interval for a given + number of samples n_samples and an accuracy gamma: + + .. math:: + \epsilon = \dfrac{1}{\sqrt{2N}}\log\left(\dfrac{2}{\gamma} \ + \right) + + Parameters + ---------- + n_samples : int + number of samples + gamma : float + accuracy + + Returns + ---------- + length of the confidence interval : float + """ + return np.sqrt(1 / (2 * n_samples) * np.log(2 / gamma)) + + @staticmethod + def confidence_intervals(a_min: float, a_max: float, ri: int): + r""" + Computes the confidence intervals + + Parameters + ---------- + a_min : float + minimum amplitude measured in the round + a_max : float + maximum amplitude measured in the round + ri : int + number of quadrants passed to get the current angle + + """ + if ri % 2 == 0: + gamma_min = np.arcsin(np.sqrt(a_min)) + gamma_max = np.arcsin(np.sqrt(a_max)) + else: + gamma_min = -np.arcsin(np.sqrt(a_max)) + 0.5 * np.pi + gamma_max = -np.arcsin(np.sqrt(a_min)) + 0.5 * np.pi + return [gamma_min, gamma_max] + + + def miqae(self, epsilon: float = 0.01, shots: int = 100, alpha: float = 0.05): + """ + This function implements Algorithm 1 from the IQAE paper. The result + is an estimation of the desired probability with precision at least + epsilon and accuracy at least alpha. + + Parameters + ---------- + epsilon : float + precision + alpha : float + accuracy + shots : int + number of measurements on each iteration + + Returns + ---------- + a_l : float + lower bound for the probability to be estimated + a_u : float + upper bound for the probability to be estimated + + """ + + self.circuit_statistics = {} + ##################################################### + i = 0 + k = int(0) + [theta_l, theta_u] = [0.0, np.pi / 2] + + # Kmax in paper + bigk_max = np.pi / (4.0 * epsilon) + + while theta_u - theta_l > 2 * epsilon: + n_ = 0 + i = i + 1 + k_i = k + # Ki in paper + big_k_i = 2 * k_i + 1 + # alpha_i in paper + alpha_i = 2.0 * alpha * big_k_i / (3.0 * bigk_max) + cte_ = (np.sin(np.pi / 21.0) * np.sin(8.0 * np.pi /21.0)) ** 2 + n_i_max = int(np.floor(2.0 * np.log(2.0 / alpha_i) / cte_)) + #number of quadrants passed + ri = np.floor(big_k_i * theta_l / (0.5 * np.pi)) + + #print("alpha_i: ", alpha_i, "n_i_max: ", n_i_max) + while k_i == k: + ##################################################### + shots_ = min(shots, n_i_max - n_) + routine = self.quantum_step(k_i) + start = time.time() + results, circuit, _, _ = get_results( + routine, + linalg_qpu=self.linalg_qpu, + shots=shots_, + qubits=self.index + ) + end = time.time() + #print("k: ", k_i, "shots: ", shots_) + self.quantum_times.append(end-start) + + if k_i not in self.schedule: + self.schedule.update({k_i:shots_}) + else: + self.schedule.update({k_i:self.schedule[k_i] + shots_}) + + n_ = n_ + shots_ + a_ = measure_state_probability(results, self.target) + epsilon_a = mIQAE.chebysev_bound(n_, alpha_i) + #print("n_: ", n_, "alpha_i: ", alpha_i, "epsilon_a: ", epsilon_a) + a_max = np.minimum(a_ + epsilon_a, 1.0) + a_min = np.maximum(a_ - epsilon_a, 0.0) + + gamma_min, gamma_max = mIQAE.confidence_intervals(a_min, a_max, ri) + + theta_l = (0.5 * np.pi * ri + gamma_min) / big_k_i + theta_u = (0.5 * np.pi * ri + gamma_max) / big_k_i + + k = mIQAE.find_next_k(k_i, theta_l, theta_u) + + [a_l, a_u] = [np.sin(theta_l) ** 2, np.sin(theta_u) ** 2] + return [a_l, a_u] + + def quantum_step(self, k): + r""" + Create the quantum routine needed for the iqae step + + Parameters + ---------- + k : int + number of Grover operator applications + + Returns + ---------- + routine : qlm routine + qlm routine for the iqae step + """ + + routine = qlm.QRoutine() + wires = routine.new_wires(self.oracle.arity) + routine.apply(self.oracle, wires) + for j in range(k): + routine.apply(self._grover_oracle, wires) + return routine + + def run(self): + r""" + run method for the class. + + Returns + ---------- + + self.ae : + amplitude estimation parameter + + """ + start = time.time() + [self.ae_l, self.ae_u] = self.miqae( + epsilon=self.epsilon, shots=self.shots, alpha=self.alpha + ) + self.theta_l = np.arcsin(np.sqrt(self.ae_l)) + self.theta_u = np.arcsin(np.sqrt(self.ae_u)) + self.theta = (self.theta_u + self.theta_l) / 2.0 + self.ae = (self.ae_u + self.ae_l) / 2.0 + end = time.time() + self.run_time = end - start + self.schedule_pdf = pd.DataFrame.from_dict( + self.schedule, + columns=['shots'], + orient='index' + ) + self.schedule_pdf.reset_index(inplace=True) + self.schedule_pdf.rename(columns={'index': 'm_k'}, inplace=True) + self.oracle_calls = np.sum( + self.schedule_pdf['shots'] * (2 * self.schedule_pdf['m_k'] + 1)) + self.max_oracle_depth = np.max(2 * self.schedule_pdf['m_k']+ 1) + self.quantum_time = sum(self.quantum_times) + return self.ae diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/modified_real_quantum_ae.py b/tnbs/BTC_02_AE/QQuantLib/AE/modified_real_quantum_ae.py new file mode 100644 index 0000000..740bfea --- /dev/null +++ b/tnbs/BTC_02_AE/QQuantLib/AE/modified_real_quantum_ae.py @@ -0,0 +1,563 @@ +""" +This module contains necessary functions and classes to implement +Real Quantum Amplitude Estimation based on the paper: + + Manzano, A., Musso, D., Leitao, A. et al. + Real Quantum Amplitude Estimation + Preprint + +Author: Gonzalo Ferro Costas & Alberto Manzano Herrero + +""" + +import time +#from copy import deepcopy +import numpy as np +import pandas as pd +import qat.lang.AQASM as qlm +from QQuantLib.AA.amplitude_amplification import grover +from QQuantLib.utils.data_extracting import get_results +from QQuantLib.utils.utils import measure_state_probability, bitfield_to_int, check_list_type, mask + + +class mRQAE: + """ + Class for Real Quantum Amplitude Estimation (RQAE) + algorithm + + Parameters + ---------- + oracle: QLM gate + QLM gate with the Oracle for implementing the + Grover operator + target : list of ints + python list with the target for the amplitude estimation + index : list of ints + qubits which mark the register to do the amplitude + estimation + + kwars : dictionary + dictionary that allows the configuration of the IQAE algorithm: \\ + Implemented keys: + + qpu : QLM solver + solver for simulating the resulting circuits + q : int + amplification ratio + epsilon : int + precision + gamma : float + accuracy + mcz_qlm : bool + for using or not QLM implementation of the multi controlled Z + gate + """ + + def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): + """ + + Method for initializing the class + + """ + ########################################### + # Setting attributes + self._oracle = oracle + self._target = check_list_type(target, int) + self._index = check_list_type(index, int) + + # Set the QPU to use + self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU + if self.linalg_qpu is None: + raise ValueError("Not QPU was provide. Please provide it!") + + self.epsilon = kwargs.get("epsilon", 0.01) + self.gamma = kwargs.get("gamma", 0.05) + # Amplification Ratio: q in the paper + self.ratio = kwargs.get("q", 2) + self.mcz_qlm = kwargs.get("mcz_qlm", True) + self.save_circuits = kwargs.get("save_circuits", False) + + # Creating the grover operator + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + self.ae_l = None + self.ae_u = None + self.ae = None + self.circuit_statistics = None + self.time_pdf = None + self.run_time = None + self.schedule = {} + self.oracle_calls = None + self.max_oracle_depth = None + self.schedule_pdf = None + self.quantum_times = [] + self.quantum_time = None + self.circuit_dict = {} + self.info = None + + ##################################################################### + @property + def oracle(self): + """ + creating oracle property + """ + return self._oracle + + @oracle.setter + def oracle(self, value): + """ + setter of the oracle property + """ + self._oracle = value + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def target(self): + """ + creating target property + """ + return self._target + + @target.setter + def target(self, value): + """ + setter of the target property + """ + self._target = check_list_type(value, int) + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def index(self): + """ + creating index property + """ + return self._index + + @index.setter + def index(self, value): + """ + setter of the index property + """ + self._index = check_list_type(value, int) + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def shifted_oracle(self): + """ + creating shifted_oracle property + """ + return self._shifted_oracle + + @shifted_oracle.setter + def shifted_oracle(self, shift): + """ + setter of the shifted_oracle property + + Parameters + ---------- + shift : float + shift for the oracle + """ + self._shifted_oracle = qlm.QRoutine() + wires = self._shifted_oracle.new_wires(self.oracle.arity + 1) + self._shifted_oracle.apply(qlm.H, wires[-1]) + self._shifted_oracle.apply( + qlm.RY(2 * np.arccos(shift)).ctrl(), wires[-1], wires[0] + ) + self._shifted_oracle.apply( + mask( + self.oracle.arity, + 2**self.oracle.arity - 1 - bitfield_to_int(self.target), + ).ctrl(), + wires[-1], + wires[: self.oracle.arity], + ) + self._shifted_oracle.apply(qlm.X, wires[-1]) + self._shifted_oracle.apply( + self._oracle.ctrl(), wires[-1], wires[: self._oracle.arity] + ) + self._shifted_oracle.apply(qlm.X, wires[-1]) + self._shifted_oracle.apply(qlm.H, wires[-1]) + + ##################################################################### + + def first_step(self, shift: float, shots: int, gamma: float): + """ + This function implements the first step of the RQAE paper. The result + is a first estimation of the desired amplitude. + + Parameters + ---------- + shift : float + shift for the first iteration + shots : int + number of measurements + gamma : float + accuracy + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + + self.shifted_oracle = 2 * shift + start = time.time() + results, circuit, _, _ = get_results( + self._shifted_oracle, self.linalg_qpu, shots=shots + ) + if self.save_circuits: + self.circuit_dict.update({"first_step": self._shifted_oracle}) + + end = time.time() + self.quantum_times.append(end-start) + start = time.time() + #probability_sum = results["Probability"].iloc[ + # bitfield_to_int([0] + list(self.target)) + #] + probability_sum = measure_state_probability( + results, [0] + list(self.target) + ) + + #probability_diff = results["Probability"].iloc[ + # bitfield_to_int([1] + list(self.target)) + #] + probability_diff = measure_state_probability( + results, [1] + list(self.target) + ) + epsilon_probability = mRQAE.chebysev_bound(shots, gamma) + + amplitude_max = np.minimum( + (probability_sum - probability_diff) / (4 * shift) + + epsilon_probability / (2 * np.abs(shift)), + 1.0, + ) + amplitude_min = np.maximum( + (probability_sum - probability_diff) / (4 * shift) + - epsilon_probability / (2 * np.abs(shift)), + -1.0, + ) + end = time.time() + first_step_time = end - start + + return [amplitude_min, amplitude_max] + + def run_step(self, shift: float, shots: int, gamma: float, k: int): + """ + This function implements a step of the RQAE paper. The result + is a refined estimation of the desired amplitude. + + Parameters + ---------- + shift : float + shift for the first iteration + shots : int + number of measurements + gamma : float + accuracy + k : int + number of amplifications + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + #print(shift) + self.shifted_oracle = 2 * shift + + grover_oracle = grover( + self.shifted_oracle, + [0] + list(self.target), + np.arange(len(self.index) + 1), + mcz_qlm=self.mcz_qlm, + ) + + routine = qlm.QRoutine() + wires = routine.new_wires(self.shifted_oracle.arity) + routine.apply(self.shifted_oracle, wires) + for i in range(k): + routine.apply(grover_oracle, wires) + start = time.time() + results, circuit, _, _ = get_results(routine, self.linalg_qpu, shots=shots) + if self.save_circuits: + self.circuit_dict.update( + {"step_{}".format(k): routine} + ) + end = time.time() + self.quantum_times.append(end-start) + #probability_sum = results["Probability"].iloc[ + # bitfield_to_int([0] + list(self.target)) + #] + probability_sum = measure_state_probability( + results, [0] + list(self.target) + ) + + epsilon_probability = mRQAE.chebysev_bound(shots, gamma) + probability_max = min(probability_sum + epsilon_probability, 1) + probability_min = max(probability_sum - epsilon_probability, 0) + angle_max = np.arcsin(np.sqrt(probability_max)) / (2 * k + 1) + angle_min = np.arcsin(np.sqrt(probability_min)) / (2 * k + 1) + amplitude_max = np.sin(angle_max) - shift + amplitude_min = np.sin(angle_min) - shift + first_step_time = end - start + + return [amplitude_min, amplitude_max] + + @staticmethod + def compute_info( + ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05, **kwargs + ): + """ + This function computes theoretical values of the IQAE algorithm. + + Parameters + ---------- + ratio: float + amplification ratio/policy + epsilon : float + precision + gamma : float + accuracy + + Return + ------ + info : dict + python dictionary with the computed information + + """ + epsilon = 0.5 * epsilon + # Maximum amplification + epsilon_infinity = 0.5 * np.sin(np.pi / (4 * ratio)) ** 2 + k_max = int( + np.ceil( + np.arcsin(np.sqrt(2 * epsilon_infinity)) + / np.arcsin(2 * epsilon) + - 0.5 + ) + ) + bigk_max = 2 * k_max + 1 + # Maximum number of grover calls + epsilon_p = 0.5 * np.sin(np.pi / (4 * (ratio + 2))) ** 2 + coef_0 = bigk_max * (ratio / (ratio - 1)) / (4 * epsilon_p ** 2) + inside_log = ratio ** (ratio / (ratio - 1)) / (gamma * (ratio - 1)) + n_grover = coef_0 * np.log(4 * np.sqrt(np.e) * inside_log) + # Shots for first iteration + epsilon_0 = 0.5 * np.sin(np.pi/(2 * (ratio + 2))) + gamma_0 = 0.5 * gamma * (ratio - 1) / (ratio * (2 * k_max + 1)) + shots_first_step = np.ceil(np.log(2 / gamma_0) / (2 * epsilon_0 ** 2)) + # Maximum number of oracle calls + n_oracle = 2 * n_grover + shots_first_step + # Maximum number of iterations + big_t = np.log( + 2.0 * ratio * ratio + * np.arcsin(np.sqrt(2.0 * epsilon_infinity)) + / np.arcsin(2.0 * epsilon) + ) / np.log(ratio) + + info = { + "epsilon_infinity": epsilon_infinity, + "epsilon_p": epsilon_p, "epsilon_0": epsilon_0, + "k_max": k_max, "n_0": shots_first_step, + "n_grover": n_grover, "n_oracle": n_oracle, + "big_t": big_t + } + + return info + @staticmethod + def display_information( + ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05, **kwargs + ): + """ + This function displays information of the properties of the + method for a given set of parameters + + Parameters + ---------- + ratio: float + amplification ratio/policy + epsilon : float + precision + gamma : float + accuracy + + """ + + + info_dict = mRQAE.compute_info( + ratio = ratio, epsilon = epsilon, gamma=gamma) + print("-------------------------------------------------------------") + print("Maximum number of amplifications: ", info_dict["k_max"]) + print("Maximum number of rounds: ", info_dict["big_t"]) + print("Maximum number of Grover operator calls: ", info_dict["n_grover"]) + print("Maximum number of Oracle operator calls: ", info_dict["n_oracle"]) + print("epsilon infinity: ", info_dict["epsilon_infinity"]) + print("Maximum number of shots for first step: ", info_dict["n_0"]) + print("-------------------------------------------------------------") + + @staticmethod + def chebysev_bound(n_samples: int, gamma: float): + """ + Computes the length of the confidence interval for a given number + of samples n_samples and an accuracy gamma. + + Parameters + ---------- + n_samples : int + number of samples + gamma : float + accuracy + + Returns + ---------- + length of the confidence interval + """ + return np.sqrt(1 / (2 * n_samples) * np.log(2 / gamma)) + + def mrqae(self, ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05): + """ + This function implements the first step of the RQAE paper. The + result is an estimation of the desired amplitude with precision + epsilon and accuracy gamma. + + Parameters + ---------- + ratio : int + amplification ratio + epsilon : int + precision + gamma : float + accuracy + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + ###################################### + + epsilon = 0.5 * epsilon + + # General parameters + epsilon_infinity = 0.5 * np.sin(np.pi / (4 * ratio)) ** 2 + k_max = int( + np.ceil( + np.arcsin(np.sqrt(2 * epsilon_infinity)) + / np.arcsin(2 * epsilon) + - 0.5 + ) + ) + + ############### First Step ####################### + # First step shift + shift_0 = 0.5 + # Bounded error first step + epsilon_0 = np.abs(shift_0) * np.sin(np.pi / (2 * (ratio + 2))) + if epsilon_0 < 2 * epsilon * np.abs(shift_0): + epsilon_0 = 2 * epsilon * np.abs(shift_0) + gamma_0 = gamma + else: + # Maximum circuit depth + # gamma for first step + gamma_0 = 0.5 * gamma * (ratio - 1) / (ratio * (2 * k_max + 1)) + # shots for first step + n_0 = int(np.ceil(np.log(2.0 / gamma_0) / (2 * epsilon_0 ** 2))) + # first step execution + [amplitude_min, amplitude_max] = self.first_step( + shift=shift_0, shots=n_0, gamma=gamma_0 + ) + self.schedule.update({0 : n_0}) + epsilon_amplitude = (amplitude_max - amplitude_min) / 2 + # print("first step. Shift ", shift_0 , "shots: ", n_0, "gamma_0: ", gamma_0) + + ############### Consecutive Steps ####################### + + while epsilon_amplitude > epsilon: + # amplification of the step + k = int(np.floor(np.pi / (4 * np.arcsin(2 * epsilon_amplitude)) - 0.5)) + # maximum k for step + k = min(k, k_max) + # new epsilon for step + epsilon_p = 0.5 * np.sin( + np.pi / (4 * (ratio + 2/(2 * k + 1))) + ) ** 2 + # shift of the step + shift = -amplitude_min + if shift > 0: + shift = min(shift, 0.5) + if shift < 0: + shift = max(shift, -0.5) + + # gamma of the step + gamma_i = 0.5 * gamma * (ratio - 1) * (2 * k + 1) \ + / (ratio * (2 * k_max + 1)) + + # number of shots of the step + n_i = int(np.ceil(np.log(2.0 / gamma_i) / (2 * epsilon_p ** 2))) + # Execute step + [amplitude_min, amplitude_max] = self.run_step( + shift=shift, shots=n_i, gamma=gamma_i, k=k + ) + # Added the shots for the k + if k not in self.schedule: + self.schedule.update({k:n_i}) + else: + # If k exists sum the shots with the before number of shots + self.schedule.update({k:self.schedule[k] + n_i}) + # print("Step k: ", k, "Shift ", shift , "shots: ", n_i, "gamma_i: ", gamma_i) + epsilon_amplitude = (amplitude_max - amplitude_min) / 2 + + return [2 * amplitude_min, 2 * amplitude_max] + + def run(self): + r""" + run method for the class. + + Returns + ---------- + + self.ae : + amplitude estimation parameter + + """ + start = time.time() + [self.ae_l, self.ae_u] = self.mrqae( + ratio=self.ratio, epsilon=self.epsilon, gamma=self.gamma + ) + # Here we write the bounds of the method + self.info = mRQAE.compute_info( + ratio=self.ratio, epsilon=self.epsilon, gamma=self.gamma + ) + self.ae = (self.ae_u + self.ae_l) / 2.0 + end = time.time() + self.run_time = end - start + self.schedule_pdf = pd.DataFrame.from_dict( + self.schedule, + columns=['shots'], + orient='index' + ) + self.schedule_pdf.reset_index(inplace=True) + self.schedule_pdf.rename(columns={'index': 'm_k'}, inplace=True) + self.oracle_calls = np.sum( + self.schedule_pdf['shots'] * (2 * self.schedule_pdf['m_k'] + 1)) + self.max_oracle_depth = np.max(2 * self.schedule_pdf['m_k']+ 1) + self.quantum_time = sum(self.quantum_times) + return self.ae diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/montecarlo_ae.py b/tnbs/BTC_02_AE/QQuantLib/AE/montecarlo_ae.py index e124f5f..e315d4f 100644 --- a/tnbs/BTC_02_AE/QQuantLib/AE/montecarlo_ae.py +++ b/tnbs/BTC_02_AE/QQuantLib/AE/montecarlo_ae.py @@ -8,20 +8,13 @@ """ -import sys -import os -import re import time +#from copy import deepcopy import numpy as np import pandas as pd import qat.lang.AQASM as qlm - - - -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.utils.data_extracting import get_results -from QQuantLib.utils.utils import check_list_type, \ - measure_state_probability +from QQuantLib.utils.utils import check_list_type, measure_state_probability class MCAE: @@ -64,9 +57,9 @@ def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): # Set the QPU to use self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") self.shots = int(kwargs.get("shots", 100)) self.mcz_qlm = kwargs.get("mcz_qlm", True) @@ -155,6 +148,8 @@ def run(self): end = time.time() self.quantum_times.append(end-start) self.ae = measure_state_probability(results, self.target) + self.ae_l = max(0.0, self.ae - 2.0 / np.sqrt(float(self.shots))) + self.ae_u = min(1.0, self.ae + 2.0 / np.sqrt(float(self.shots))) self.run_time = end - start self.schedule_pdf = pd.DataFrame( [[0, self.shots]], diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/real_quantum_ae.py b/tnbs/BTC_02_AE/QQuantLib/AE/real_quantum_ae.py index 5557e0e..c15f164 100644 --- a/tnbs/BTC_02_AE/QQuantLib/AE/real_quantum_ae.py +++ b/tnbs/BTC_02_AE/QQuantLib/AE/real_quantum_ae.py @@ -10,20 +10,14 @@ """ -import sys -import os -import re import time +#from copy import deepcopy import numpy as np import pandas as pd import qat.lang.AQASM as qlm - - -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.AA.amplitude_amplification import grover from QQuantLib.utils.data_extracting import get_results -from QQuantLib.utils.utils import measure_state_probability, \ - bitfield_to_int, check_list_type, mask +from QQuantLib.utils.utils import measure_state_probability, bitfield_to_int, check_list_type, mask class RQAE: @@ -73,9 +67,9 @@ def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): # Set the QPU to use self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") self.epsilon = kwargs.get("epsilon", 0.01) self.gamma = kwargs.get("gamma", 0.05) @@ -93,6 +87,7 @@ def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): self.ae_u = None self.ae = None self.circuit_statistics = None + self.circuit_statistics = {} self.time_pdf = None self.run_time = None self.schedule = {} @@ -102,6 +97,7 @@ def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): self.quantum_times = [] self.quantum_time = None self.circuit_dict = {} + self.info = None ##################################################################### @property @@ -229,10 +225,6 @@ def first_step(self, shift: float, shots: int, gamma: float): end = time.time() self.quantum_times.append(end-start) start = time.time() - step_circuit_stats = circuit.statistics() - step_circuit_stats.update({"n_shots": shots}) - self.circuit_statistics.update({0: step_circuit_stats}) - self.schedule.update({0 : shots}) #probability_sum = results["Probability"].iloc[ # bitfield_to_int([0] + list(self.target)) @@ -262,7 +254,7 @@ def first_step(self, shift: float, shots: int, gamma: float): end = time.time() first_step_time = end - start - return [amplitude_min, amplitude_max] + return [amplitude_min, amplitude_max], circuit def run_step(self, shift: float, shots: int, gamma: float, k: int): """ @@ -311,10 +303,6 @@ def run_step(self, shift: float, shots: int, gamma: float, k: int): ) end = time.time() self.quantum_times.append(end-start) - step_circuit_stats = circuit.statistics() - step_circuit_stats.update({"n_shots": shots}) - self.circuit_statistics.update({k: step_circuit_stats}) - self.schedule.update({k : shots}) #probability_sum = results["Probability"].iloc[ # bitfield_to_int([0] + list(self.target)) #] @@ -331,15 +319,14 @@ def run_step(self, shift: float, shots: int, gamma: float, k: int): amplitude_min = np.sin(angle_min) - shift first_step_time = end - start - return [amplitude_min, amplitude_max] + return [amplitude_min, amplitude_max], circuit @staticmethod - def display_information( - ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05 + def compute_info( + ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05, **kwargs ): """ - This function displays information of the properties of the - method for a given set of parameters + This function computes theoretical values of the IQAE algorithm. Parameters ---------- @@ -350,34 +337,79 @@ def display_information( gamma : float accuracy + Return + ------ + info : dict + python dictionary with the computed information + """ - theoretical_epsilon = 0.5 * np.sin(np.pi / (2 * (ratio + 2))) ** 2 + epsilon = 0.5 * epsilon + # Bounded for the error at each step + theoretical_epsilon = 0.5 * np.sin(np.pi / (4.0 * (ratio + 2))) ** 2 + # Maximum amplification k_max = int( np.ceil( np.arcsin(np.sqrt(2 * theoretical_epsilon)) / np.arcsin(2 * epsilon) - * 0.5 - 0.5 ) ) bigk_max = 2 * k_max + 1 + # Maximum number of iterations big_t = np.log( - ratio + 2.0 + * ratio * ratio * (np.arcsin(np.sqrt(2 * theoretical_epsilon))) / (np.arcsin(2 * epsilon)) ) / np.log(ratio) + # Maximum probability failure at each step gamma_i = gamma / big_t # This is shots for each iteration: Ni in the paper n_i = int( np.ceil(1 / (2 * theoretical_epsilon**2) * np.log(2 * big_t / gamma)) ) - n_oracle = int(n_i / 2 * bigk_max * (1 + ratio / (ratio - 1))) + # Total number of Grover operator calls + n_grover = int(n_i / 2 * bigk_max * (1 + ratio / (ratio - 1))) + # This is the number of calls to the oracle operator (A) + n_oracle = 2 * n_grover + n_i + + info = { + "theoretical_epsilon": theoretical_epsilon, "k_max": k_max, + "big_t": big_t, "gamma_i": gamma_i, "n_i": n_i, + "n_grover": n_grover, "n_oracle": n_oracle, + } + + return info + + @staticmethod + def display_information( + ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05, **kwargs + ): + """ + This function displays information of the properties of the + method for a given set of parameters + + Parameters + ---------- + ratio: float + amplification ratio/policy + epsilon : float + precision + gamma : float + accuracy + + """ + + info_dict = RQAE.compute_info( + ratio = ratio, epsilon = epsilon, gamma=gamma) + print("-------------------------------------------------------------") - print("Maximum number of amplifications: ", k_max) - print("Maximum number of rounds: ", int(big_t)) - print("Number of shots per round: ", n_i) - print("Maximum number of calls to the oracle: ", n_oracle) + print("Maximum number of amplifications: ", info_dict["k_max"]) + print("Maximum number of rounds: ", info_dict["big_t"]) + print("Number of shots per round: ", info_dict["n_i"]) + print("Maximum number of Grover operator calls: ", info_dict["n_grover"]) + print("Maximum number of Oracle operator calls: ", info_dict["n_oracle"]) print("-------------------------------------------------------------") @staticmethod @@ -426,20 +458,18 @@ def rqae(self, ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05): epsilon = 0.5 * epsilon # Always need to clean the circuit statistics property - self.circuit_statistics = {} - # time_list = [] - theoretical_epsilon = 0.5 * np.sin(np.pi / (2 * (ratio + 2))) ** 2 + theoretical_epsilon = 0.5 * np.sin(np.pi / (4.0 * (ratio + 2))) ** 2 k_max = int( np.ceil( np.arcsin(np.sqrt(2 * theoretical_epsilon)) / np.arcsin(2 * epsilon) - * 0.5 - 0.5 ) ) bigk_max = 2 * k_max + 1 big_t = np.log( - ratio + 2.0 + * ratio * ratio * (np.arcsin(np.sqrt(2 * theoretical_epsilon))) / (np.arcsin(2 * epsilon)) @@ -451,24 +481,38 @@ def rqae(self, ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05): ) epsilon_probability = np.sqrt(1 / (2 * n_i) * np.log(2 / gamma_i)) shift = theoretical_epsilon / np.sin(np.pi / (2 * (ratio + 2))) - #print("First Step: {}".format(shift)) + # print("first step. Shift ", shift , "shots: ", n_i, "gamma_0: ", gamma_i) + ##################################### # First step - [amplitude_min, amplitude_max] = self.first_step( + [amplitude_min, amplitude_max], _ = self.first_step( shift=shift, shots=n_i, gamma=gamma_i ) epsilon_amplitude = (amplitude_max - amplitude_min) / 2 + # Added step to schedule: m_k, shots + self.schedule.update({0 : n_i}) + # time_list.append(time_pdf) # Consecutive steps while epsilon_amplitude > epsilon: k = int(np.floor(np.pi / (4 * np.arcsin(2 * epsilon_amplitude)) - 0.5)) k = min(k, k_max) shift = -amplitude_min - shift = min(shift, 0.5) - #print("While Step: {}".format(shift)) - [amplitude_min, amplitude_max] = self.run_step( + if shift > 0: + shift = min(shift, 0.5) + if shift < 0: + shift = max(shift, -0.5) + # print("Step k: ", k, "Shift ", shift , "shots: ", n_i, "gamma_0: ", gamma_i) + [amplitude_min, amplitude_max], _ = self.run_step( shift=shift, shots=n_i, gamma=gamma_i, k=k ) + # Added the shots for the k + if k not in self.schedule: + self.schedule.update({k:n_i}) + else: + # If k exists sum the shots with the before number of shots + self.schedule.update({k:self.schedule[k] + n_i}) + # time_list.append(time_pdf) epsilon_amplitude = (amplitude_max - amplitude_min) / 2 @@ -489,6 +533,10 @@ def run(self): [self.ae_l, self.ae_u] = self.rqae( ratio=self.ratio, epsilon=self.epsilon, gamma=self.gamma ) + # Here we write the bounds of the method + self.info = RQAE.compute_info( + ratio=self.ratio, epsilon=self.epsilon, gamma=self.gamma + ) self.ae = (self.ae_u + self.ae_l) / 2.0 end = time.time() self.run_time = end - start diff --git a/tnbs/BTC_02_AE/QQuantLib/AE/shots_real_quantum_ae.py b/tnbs/BTC_02_AE/QQuantLib/AE/shots_real_quantum_ae.py new file mode 100644 index 0000000..e20a3da --- /dev/null +++ b/tnbs/BTC_02_AE/QQuantLib/AE/shots_real_quantum_ae.py @@ -0,0 +1,569 @@ +""" +This module contains necessary functions and classes to implemeneal Quantum Amplitude Estimation based on the paper: + + Manzano, A., Musso, D., Leitao, A. et al. + Real Quantum Amplitude Estimation + Preprint + +Author: Gonzalo Ferro Costas & Alberto Manzano Herrero + +""" + +import time +#from copy import deepcopy +import numpy as np +import pandas as pd +import qat.lang.AQASM as qlm +from QQuantLib.AA.amplitude_amplification import grover +from QQuantLib.utils.data_extracting import get_results +from QQuantLib.utils.utils import measure_state_probability, bitfield_to_int, check_list_type, mask + + +class sRQAE: + """ + Class for Real Quantum Amplitude Estimation (RQAE) + algorithm + + Parameters + ---------- + oracle: QLM gate + QLM gate with the Oracle for implementing the + Grover operator + target : list of ints + python list with the target for the amplitude estimation + index : list of ints + qubits which mark the register to do the amplitude + estimation + + kwars : dictionary + dictionary that allows the configuration of the IQAE algorithm: \\ + Implemented keys: + + qpu : QLM solver + solver for simulating the resulting circuits + q : int + amplification ratio + epsilon : int + precision + gamma : float + accuracy + mcz_qlm : bool + for using or not QLM implementation of the multi controlled Z + gate + """ + + def __init__(self, oracle: qlm.QRoutine, target: list, index: list, **kwargs): + """ + + Method for initializing the class + + """ + ########################################### + # Setting attributes + self._oracle = oracle + self._target = check_list_type(target, int) + self._index = check_list_type(index, int) + + # Set the QPU to use + self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU + if self.linalg_qpu is None: + raise ValueError("Not QPU was provide. Please provide it!") + + self.epsilon = kwargs.get("epsilon", 0.01) + self.gamma = kwargs.get("gamma", 0.05) + # Amplification Ratio: q in the paper + self.ratio = kwargs.get("q", 3) + self.shots = int(kwargs.get("shots", 100)) + self.mcz_qlm = kwargs.get("mcz_qlm", True) + self.save_circuits = kwargs.get("save_circuits", False) + + # Creating the grover operator + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + self.ae_l = None + self.ae_u = None + self.ae = None + self.circuit_statistics = None + self.circuit_statistics = {} + self.time_pdf = None + self.run_time = None + self.schedule = {} + self.oracle_calls = None + self.max_oracle_depth = None + self.schedule_pdf = None + self.quantum_times = [] + self.quantum_time = None + self.info = None + self.circuit_dict = {} + + ##################################################################### + @property + def oracle(self): + """ + creating oracle property + """ + return self._oracle + + @oracle.setter + def oracle(self, value): + """ + setter of the oracle property + """ + self._oracle = value + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def target(self): + """ + creating target property + """ + return self._target + + @target.setter + def target(self, value): + """ + setter of the target property + """ + self._target = check_list_type(value, int) + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def index(self): + """ + creating index property + """ + return self._index + + @index.setter + def index(self, value): + """ + setter of the index property + """ + self._index = check_list_type(value, int) + self._grover_oracle = grover( + self._oracle, self.target, self.index, mcz_qlm=self.mcz_qlm + ) + + @property + def shifted_oracle(self): + """ + creating shifted_oracle property + """ + return self._shifted_oracle + + @shifted_oracle.setter + def shifted_oracle(self, shift): + """ + setter of the shifted_oracle property + + Parameters + ---------- + shift : float + shift for the oracle + """ + self._shifted_oracle = qlm.QRoutine() + wires = self._shifted_oracle.new_wires(self.oracle.arity + 1) + self._shifted_oracle.apply(qlm.H, wires[-1]) + self._shifted_oracle.apply( + qlm.RY(2 * np.arccos(shift)).ctrl(), wires[-1], wires[0] + ) + self._shifted_oracle.apply( + mask( + self.oracle.arity, + 2**self.oracle.arity - 1 - bitfield_to_int(self.target), + ).ctrl(), + wires[-1], + wires[: self.oracle.arity], + ) + self._shifted_oracle.apply(qlm.X, wires[-1]) + self._shifted_oracle.apply( + self._oracle.ctrl(), wires[-1], wires[: self._oracle.arity] + ) + self._shifted_oracle.apply(qlm.X, wires[-1]) + self._shifted_oracle.apply(qlm.H, wires[-1]) + + ##################################################################### + + def first_step(self, shift: float, shots: int, gamma: float): + """ + This function implements the first step of the RQAE paper. The result + is a first estimation of the desired amplitude. + + Parameters + ---------- + shift : float + shift for the first iteration + shots : int + number of measurements + gamma : float + accuracy + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + + self.shifted_oracle = 2 * shift + start = time.time() + results, circuit, _, _ = get_results( + self._shifted_oracle, self.linalg_qpu, shots=shots + ) + if self.save_circuits: + self.circuit_dict.update({"first_step": self._shifted_oracle}) + + end = time.time() + self.quantum_times.append(end-start) + + probability_sum = measure_state_probability( + results, [0] + list(self.target) + ) + + probability_diff = measure_state_probability( + results, [1] + list(self.target) + ) + return probability_sum, probability_diff + + def run_step(self, shift: float, shots: int, gamma: float, k: int): + """ + This function implements a step of the RQAE paper. The result + is a refined estimation of the desired amplitude. + + Parameters + ---------- + shift : float + shift for the first iteration + shots : int + number of measurements + gamma : float + accuracy + k : int + number of amplifications + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + #print(shift) + self.shifted_oracle = 2 * shift + + grover_oracle = grover( + self.shifted_oracle, + [0] + list(self.target), + np.arange(len(self.index) + 1), + mcz_qlm=self.mcz_qlm, + ) + + routine = qlm.QRoutine() + wires = routine.new_wires(self.shifted_oracle.arity) + routine.apply(self.shifted_oracle, wires) + for i in range(k): + routine.apply(grover_oracle, wires) + start = time.time() + results, circuit, _, _ = get_results(routine, self.linalg_qpu, shots=shots) + if self.save_circuits: + self.circuit_dict.update( + {"step_{}".format(k): routine} + ) + end = time.time() + self.quantum_times.append(end-start) + #probability_sum = results["Probability"].iloc[ + # bitfield_to_int([0] + list(self.target)) + #] + probability_sum = measure_state_probability( + results, [0] + list(self.target) + ) + return probability_sum + + + @staticmethod + def compute_info( + ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05, shots: int = 100 + ): + """ + This function computes theoretical values of the IQAE algorithm. + + Parameters + ---------- + ratio: float + amplification ratio/policy + epsilon : float + precision + gamma : float + accuracy + + Return + ------ + info : dict + python dictionary with the computed information + + """ + epsilon = 0.5 * epsilon + # Bounded for the error at each step + #theoretical_epsilon = 0.5 * np.sin(np.pi / (2 * (ratio + 2))) ** 2 + theoretical_epsilon = None + k_max = None + bigk_max = None + big_t = None + # Maximum probability failure at each step + gamma_i = None + # This is shots for each iteration: Ni in the paper + shots_max = None + # Total number of Grover operator calls + n_grover = None + # This is the number of calls to the oracle operator (A) + n_oracle = None + + info = { + "theoretical_epsilon": theoretical_epsilon, "k_max": k_max, + "big_t": big_t, "gamma_i": gamma_i, "shots_max": shots_max, + "n_grover": n_grover, "n_oracle": n_oracle, + } + + return info + + @staticmethod + def display_information( + ratio: float = 2, epsilon: float = 0.01, gamma: float = 0.05, shots: int = 100 + ): + """ + This function displays information of the properties of the + method for a given set of parameters + + Parameters + ---------- + ratio: float + amplification ratio/policy + epsilon : float + precision + gamma : float + accuracy + + """ + + info_dict = sRQAE.compute_info( + ratio = ratio, epsilon=epsilon, gamma=gamma) + + print("-------------------------------------------------------------") + print("BE AWARE: In RQAE with shots the bounds depend on the shots") + print("Here Not info is provided.") + print("Maximum number of amplifications: ", info_dict["k_max"]) + print("Maximum number of rounds: ", info_dict["big_t"]) + print("Maximum number of Grover operator calls: ", info_dict["n_grover"]) + print("Maximum number of Oracle operator calls: ", info_dict["n_oracle"]) + print("-------------------------------------------------------------") + + @staticmethod + def chebysev_bound(n_samples: int, gamma: float): + """ + Computes the length of the confidence interval for a given number + of samples n_samples and an accuracy gamma. + + Parameters + ---------- + n_samples : int + number of samples + gamma : float + accuracy + + Returns + ---------- + length of the confidence interval + """ + return np.sqrt(1 / (2 * n_samples) * np.log(2 / gamma)) + + def srqae(self, epsilon: float = 0.01, gamma: float = 0.05, user_shots: int = 100, ratio: float = 2.0): + """ + This function implements the first step of the RQAE paper. The + result is an estimation of the desired amplitude with precision + epsilon and accuracy gamma. + + Parameters + ---------- + ratio : int + amplification ratio + epsilon : int + precision + gamma : float + accuracy + shots : int + shots + + Returns + ---------- + amplitude_min : float + lower bound for the amplitude to be estimated + amplitude_max : float + upper bound for the amplitude to be estimated + + """ + ###################################### + + epsilon = 0.5 * epsilon + # time_list = [] + theoretical_epsilon = 0.5 * np.sin(np.pi / (4.0 * (ratio + 2))) ** 2 + k_max = int( + np.ceil( + #np.arcsin(np.sqrt(2 * theoretical_epsilon)) + 0.5 * np.pi + / np.arcsin(2 * epsilon) + - 0.5 + ) + ) + bigk_max = 2 * k_max + 1 + big_t = np.ceil(np.log( + ratio + * ratio + #* (np.arcsin(np.sqrt(2 * theoretical_epsilon))) + * np.pi + / (np.arcsin(2 * epsilon)) + ) / np.log(ratio)) + gamma_i = gamma / big_t + shots_max = int(np.ceil(1 / (2 * theoretical_epsilon ** 2) * np.log(2 * big_t / gamma))) + + # Compute shift for first iteration + shift = theoretical_epsilon / np.sin(np.pi / (2 * (ratio + 2))) + # Compute shots for first iteration + shots = min(shots_max, user_shots) + + ################ FIRST STEP ################### + h_plus = 0 + h_minus = 0 + n_effective = 0 + k = 0 + big_k = 2 * k + 1 + epsilon_amplitude = 0.5 + while ((big_k < ratio) and (epsilon_amplitude > epsilon)): + # First step: return probability of sum and rest + p_plus, p_minus = self.first_step( + shift=shift, shots=shots, gamma=gamma_i + ) + + # combine iterations + h_plus = h_plus + int(p_plus * shots) + h_minus = h_minus + int(p_minus * shots) + n_effective = n_effective + shots + p_plus = h_plus / n_effective + p_minus = h_minus / n_effective + epsilon_probability = sRQAE.chebysev_bound(n_effective, gamma_i) + amplitude_max = np.minimum( + (p_plus - p_minus) / (4 * shift) + + epsilon_probability / (2 * np.abs(shift)), + 0.5, + ) + amplitude_min = np.maximum( + (p_plus - p_minus) / (4 * shift) + - epsilon_probability / (2 * np.abs(shift)), + -0.5, + ) + epsilon_amplitude = (amplitude_max - amplitude_min) / 2 + self.schedule.update({k:n_effective}) + k = int(np.floor(np.pi / (4 * np.arcsin(2 * epsilon_amplitude)) - 0.5)) + big_k = 2 * k + 1 + # print( + # "first step. Shift ", shift , "shots: ", shots, "n_effective: ", n_effective, + # "amplitude_max", amplitude_max, "amplitude_min", amplitude_min + # ) + shots = min(shots_max - n_effective, user_shots) + + ################ NEXT STEPS ################### + + k_old = 0 + big_k_old = 2 * k_old + 1 + h_k = 0 + n_effective = 0 + shots = min(shots_max, user_shots) + while epsilon_amplitude > epsilon: + k = int(np.floor(np.pi / (4 * np.arcsin(2 * epsilon_amplitude)) - 0.5)) + k = min(k, k_max) + big_k = 2 *k + 1 + + if big_k < big_k_old * ratio: + k = k_old + #print("DENTRO!") + shots = min(shots_max - n_effective, user_shots) + else: + h_k = 0 + n_effective = 0 + shots = min(shots_max, user_shots) + shift = -amplitude_min + if shift > 0: + shift = min(shift, 0.5) + if shift < 0: + shift = max(shift, -0.5) + p_sum = self.run_step( + shift=shift, shots=shots, gamma=gamma_i, k=k + ) + h_k = h_k + int(p_sum * shots) + n_effective = n_effective + shots + p_sum = h_k / n_effective + self.schedule.update({k:n_effective}) + + epsilon_probability = sRQAE.chebysev_bound(n_effective, gamma_i) + probability_max = min(p_sum + epsilon_probability, 1) + probability_min = max(p_sum - epsilon_probability, 0) + angle_max = np.arcsin(np.sqrt(probability_max)) / (2 * k + 1) + angle_min = np.arcsin(np.sqrt(probability_min)) / (2 * k + 1) + amplitude_max = np.sin(angle_max) - shift + amplitude_min = np.sin(angle_min) - shift + # print( + # "Step k: ", k, "Shift ", shift , "shots: ", shots, "n_effective: ", n_effective, + # "amplitude_max", amplitude_max, "amplitude_min", amplitude_min + # ) + epsilon_amplitude = (amplitude_max - amplitude_min) / 2 + k_old = k + big_k_old = 2 * k_old + 1 + + return [2 * amplitude_min, 2 * amplitude_max] + + def run(self): + r""" + run method for the class. + + Returns + ---------- + + self.ae : + amplitude estimation parameter + + """ + start = time.time() + # print("RUN") + # print(self.epsilon, self.gamma, self.shots, self.ratio) + # print("RUN") + [self.ae_l, self.ae_u] = self.srqae( + epsilon=self.epsilon, gamma=self.gamma, user_shots=self.shots, ratio=self.ratio + ) + self.info = self.compute_info( + epsilon=self.epsilon, gamma=self.gamma, ratio=self.ratio, shots=self.shots + ) + + self.ae = (self.ae_u + self.ae_l) / 2.0 + end = time.time() + self.run_time = end - start + self.schedule_pdf = pd.DataFrame.from_dict( + self.schedule, + columns=['shots'], + orient='index' + ) + self.schedule_pdf.reset_index(inplace=True) + self.schedule_pdf.rename(columns={'index': 'm_k'}, inplace=True) + self.oracle_calls = np.sum( + self.schedule_pdf['shots'] * (2 * self.schedule_pdf['m_k'] + 1)) + self.max_oracle_depth = np.max(2 * self.schedule_pdf['m_k']+ 1) + self.quantum_time = sum(self.quantum_times) + return self.ae diff --git a/tnbs/BTC_02_AE/QQuantLib/DL/data_loading.py b/tnbs/BTC_02_AE/QQuantLib/DL/data_loading.py index 91aac28..7b9e6ee 100644 --- a/tnbs/BTC_02_AE/QQuantLib/DL/data_loading.py +++ b/tnbs/BTC_02_AE/QQuantLib/DL/data_loading.py @@ -18,14 +18,10 @@ """ -import sys -import os import time -import re import numpy as np import qat.lang.AQASM as qlm -from QQuantLib.utils.utils import mask, fwht, \ - left_conditional_probability, expmod +from QQuantLib.utils.utils import mask, fwht, left_conditional_probability, expmod # Loading uniform distribution @qlm.build_gate("UD", [int], arity=lambda x: x) @@ -202,7 +198,7 @@ def load_angles(angles: np.array, method: str = "multiplexor"): def load_array( function_array: np.array, method: str = "multiplexor", - id_name: str = str(time.time_ns()), + id_name: str = None, ): """ Creates a QLM AbstractGate for loading a normalised array into a quantum @@ -228,6 +224,9 @@ def load_array( """ number_qubits = int(np.log2(function_array.size)) + 1 + if id_name is None: + id_name = str(time.time_ns()) + @qlm.build_gate("F_{" + id_name + "}", [], arity=number_qubits) def load_array_gate(): """ @@ -245,7 +244,8 @@ def load_array_gate(): def load_probability( probability_array: np.array, method: str = "multiplexor", - id_name: str = str(time.time_ns()) + id_name: str = None, + #id_name: str = str(time.time_ns()) ): """ Creates a QLM Abstract gate for loading a given discretized probability @@ -271,6 +271,8 @@ def load_probability( Quantum Multiplexors """ number_qubits = int(np.log2(probability_array.size)) + if id_name is None: + id_name = str(time.time_ns()) @qlm.build_gate("P_{" + id_name + "}", [], arity=number_qubits) def load_probability_gate(): diff --git a/tnbs/BTC_02_AE/QQuantLib/DL/encoding_protocols.py b/tnbs/BTC_02_AE/QQuantLib/DL/encoding_protocols.py index a5e61e2..3a49600 100644 --- a/tnbs/BTC_02_AE/QQuantLib/DL/encoding_protocols.py +++ b/tnbs/BTC_02_AE/QQuantLib/DL/encoding_protocols.py @@ -5,9 +5,6 @@ """ -import sys -import os -import re import warnings import numpy as np import qat.lang.AQASM as qlm @@ -167,7 +164,9 @@ def oracle_encoding_0(self): warnings.warn('Some elements of the input array_function are negative') # Creation of function loading gate self.function_gate = dl.load_array( - np.sqrt(np.abs(self.function)), id_name="Function", method=self.multiplexor) + np.sqrt(np.abs(self.function)), + method=self.multiplexor + ) self.registers = self.oracle.new_wires(self.function_gate.arity) # Step 1 of Procedure: apply loading probability gate self.oracle.apply(self.p_gate, self.registers[: self.p_gate.arity]) @@ -223,7 +222,6 @@ def oracle_encoding_1(self): # Step 3 of Procedure: apply loading function operator for loading p(x) self.p_gate = dl.load_array( self.probability, - id_name="Probability", method=self.multiplexor ) self.oracle.apply( @@ -232,7 +230,6 @@ def oracle_encoding_1(self): # Step 5 of Procedure: apply loading function operator for loading f(x) self.function_gate = dl.load_array( self.function, - id_name="Function", method=self.multiplexor ) self.oracle.apply( @@ -280,7 +277,6 @@ def oracle_encoding_2(self): if self.probability is not None: self.p_gate = dl.load_probability( self.probability, - id_name="Probability", method=self.multiplexor ) else: @@ -289,7 +285,6 @@ def oracle_encoding_2(self): # Creation of function loading gate self.function_gate = dl.load_array( self.function, - id_name="Function", method=self.multiplexor ) self.registers = self.oracle.new_wires(self.function_gate.arity) @@ -302,10 +297,8 @@ def oracle_encoding_2(self): self.target = [0 for i in range(self.oracle.arity)] self.index = [i for i in range(self.oracle.arity)] + def run(self): - """ - Execute the computations - """ if self.encoding is None: error_string = ( "Encoding parameter MUST NOT BE None." diff --git a/tnbs/BTC_02_AE/QQuantLib/PE/classical_qpe.py b/tnbs/BTC_02_AE/QQuantLib/PE/classical_qpe.py index 8ee7066..99b3bea 100644 --- a/tnbs/BTC_02_AE/QQuantLib/PE/classical_qpe.py +++ b/tnbs/BTC_02_AE/QQuantLib/PE/classical_qpe.py @@ -15,13 +15,9 @@ """ -import sys -import os -import re import time +import numpy as np import qat.lang.AQASM as qlm - -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.utils.utils import load_qn_gate from QQuantLib.utils.data_extracting import get_results @@ -72,9 +68,9 @@ def __init__(self, **kwargs): # Set the QPU to use self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") self.shots = kwargs.get("shots", 10) self.complete = kwargs.get("complete", False) diff --git a/tnbs/BTC_02_AE/QQuantLib/PE/iterative_quantum_pe.py b/tnbs/BTC_02_AE/QQuantLib/PE/iterative_quantum_pe.py index 79a2612..b988062 100644 --- a/tnbs/BTC_02_AE/QQuantLib/PE/iterative_quantum_pe.py +++ b/tnbs/BTC_02_AE/QQuantLib/PE/iterative_quantum_pe.py @@ -14,15 +14,12 @@ """ -import sys -import os -import re import time +#from copy import deepcopy import numpy as np import pandas as pd import qat.lang.AQASM as qlm from qat.core import Result -from QQuantLib.utils.qlm_solver import get_qpu from QQuantLib.utils.data_extracting import ( create_qprogram, create_qjob, @@ -83,9 +80,9 @@ def __init__(self, **kwargs): # Set the QPU to use self.linalg_qpu = kwargs.get("qpu", None) + # Provide QPU if self.linalg_qpu is None: - print("Not QPU was provide. PyLinalg will be used") - self.linalg_qpu = get_qpu("python") + raise ValueError("Not QPU was provide. Please provide it!") self.shots = kwargs.get("shots", 10) # self.zalo = kwargs.get('zalo', False) diff --git a/tnbs/BTC_02_AE/QQuantLib/finance/ae_price_estimation.py b/tnbs/BTC_02_AE/QQuantLib/finance/ae_price_estimation.py deleted file mode 100644 index 22eb92f..0000000 --- a/tnbs/BTC_02_AE/QQuantLib/finance/ae_price_estimation.py +++ /dev/null @@ -1,162 +0,0 @@ -""" -Function for automation of different option price estimation using -Amplitude Estimation algorithms. - -The function uses the DensityProbability and the PayOff classes for -defining the option price estimation problem. Then q_solve_integral -function is used for computing the expected value integral. The function -deals with all the mandatory normalisations for returned the desired price -estimation. - -Authors: Alberto Pedro Manzano Herrero & Gonzalo Ferro Costas -""" - -import numpy as np -import pandas as pd - -from probability_class import DensityProbability -from payoff_class import PayOff -from quantum_integration import q_solve_integral - - -def ae_price_estimation(**kwargs): - """ - Configures an option price estimation problem and solving it using - AE integration techniques - - Parameters - ---------- - - kwargs : dictionary. - Dictionary for configuring the price estimation problem, the - encoding of the price estimation data into the quantum circuit - and the AE integration technique for solving it. - - Note - ____ - - The keys for the input kwargs dictionary will be the necessary keys - for configuring the DensityProbability class \\ - (see QQuantLib.finance.probability_class), the PayOff class \\ - (see QQuantLib.finance.payoff_class) and the q_solve_integral \\ - function (see QQuantLib.finance.quantum_integration). - - Returns - _______ - - pdf : Pandas DataFrame - DataFrame with the configuration of the AE problem and the solution - """ - - ae_problem = kwargs - #Building the domain - n_qbits = ae_problem.get("n_qbits", None) - x0 = ae_problem.get("x0", 1.0) - xf = ae_problem.get("xf", 3.0) - domain = np.linspace(x0, xf, 2**n_qbits) - - #Building the Probability distribution - pc = DensityProbability(**ae_problem) - p_x = pc.probability(domain) - #Normalisation of the probability distribution - p_x_normalisation = np.sum(p_x) + 1e-8 - norm_p_x = p_x / p_x_normalisation - - #Building the option payoff - po = PayOff(**ae_problem) - pay_off = po.pay_off(domain) - #Normalisation of the pay off - pay_off_normalisation = np.max(np.abs(pay_off)) + 1e-8 - norm_pay_off = pay_off / pay_off_normalisation - - #Getting the exact price of the option under BS - exact_solution = None - if po.pay_off_bs is not None: - exact_solution = po.pay_off_bs(**ae_problem) - - lista = [] - #For doing several repetitions - for i in range(ae_problem["number_of_tests"]): - #Each loop step solves a complete price estimation problem - - #Now we update the input dictionary with the probability and the - #function arrays - ae_problem.update({ - "array_function" : norm_pay_off, - "array_probability" : norm_p_x, - }) - - #EXECUTE COMPUTATION - solution, solver_object = q_solve_integral(**ae_problem) - - #For generating the output DataFrame we delete the arrays - del ae_problem["array_function"] - del ae_problem["array_probability"] - - #Undoing the normalisations - ae_expectation = solution * pay_off_normalisation * p_x_normalisation - - #Creating the output DataFrame with the complete information - - #The basis will be the input python dictionary for traceability - pdf = pd.DataFrame([ae_problem]) - #Added normalisation constants - pdf["payoff_normalisation"] = pay_off_normalisation - pdf["p_x_normalisation"] = p_x_normalisation - - #Expectation calculation using Riemann sum - pdf["riemann_expectation"] = np.sum(p_x * pay_off) - #Expectation calculation using AE integration techniques - pdf[ - [col + "_expectation" for col in ae_expectation.columns] - ] = ae_expectation - - #Option price estimation using expectation computed as Riemann sum - pdf["riemann_price_estimation"] = pdf["riemann_expectation"] * np.exp( - -pdf["risk_free_rate"] * pdf["maturity"] - ) - #Exact option price under the Black-Scholes model - pdf["exact_price"] = exact_solution - #Option price estimation using expectation computed by AE integration - pdf[[col + "_price_estimation" for col in ae_expectation.columns]] = ( - ae_expectation - * np.exp(-pdf["risk_free_rate"] * pdf["maturity"]).iloc[0] - ) - #Computing Absolute: Rieman vs AE techniques - pdf["error_riemann"] = np.abs( - pdf["ae_price_estimation"] - pdf["riemann_price_estimation"] - ) - #Computing Relative: Rieman vs AE techniques - pdf["relative_error_riemann"] = ( - pdf["error_riemann"] / pdf["riemann_price_estimation"] - ) - #Computing Absolute error: Exact BS price vs AE techniques - pdf["error_exact"] = np.abs( - pdf["ae_price_estimation"] - pdf["exact_price"]) - #Computing Relative error: Exact BS price vs AE techniques - pdf["relative_error_exact"] = pdf["error_exact"] / pdf["exact_price"] - #Other interesting staff - if solver_object is None: - #Computation Fails Encoding 0 and RQAE - pdf["schedule_pdf"] = [None] - pdf["oracle_calls"] = [None] - pdf["max_oracle_depth"] = [None] - pdf["circuit_stasts"] = [None] - pdf["run_time"] = [None] - else: - if solver_object.schedule_pdf is None: - pdf["schedule_pdf"] = [None] - else: - pdf["schedule_pdf"] = [solver_object.schedule_pdf.to_dict()] - pdf["oracle_calls"] = solver_object.oracle_calls - pdf["max_oracle_depth"] = solver_object.max_oracle_depth - pdf["circuit_stasts"] = [solver_object.solver_ae.circuit_statistics] - pdf["run_time"] = solver_object.solver_ae.run_time - - #Saving pdf - if ae_problem["save"]: - with open(ae_problem["file_name"], "a") as f_pointer: - pdf.to_csv(f_pointer, mode="a", header=f_pointer.tell() == 0, sep = ";") - lista.append(pdf) - complete_pdf = pd.concat(lista) - return complete_pdf diff --git a/tnbs/BTC_02_AE/QQuantLib/finance/classical_finance.py b/tnbs/BTC_02_AE/QQuantLib/finance/classical_finance.py deleted file mode 100644 index be3489f..0000000 --- a/tnbs/BTC_02_AE/QQuantLib/finance/classical_finance.py +++ /dev/null @@ -1,691 +0,0 @@ -""" -This module contains several auxiliary functions used in quantitative -classical finances - -Authors: Alberto Pedro Manzano Herrero & Gonzalo Ferro Costas -""" - -import numpy as np -from scipy.stats import norm - -def call_payoff(s_t: float, strike: float, **kwargs): - r"""Computes the payoff of a European call option. - - Notes - ----- - - .. math:: - C(S_T, K) = \left(S_T-K, 0\right)^+ - - Parameters - ---------- - s_t : float - price - strike : float - the strike - - Returns - ------- - payoff : float - the payoff - """ - return np.maximum(s_t - strike, 0) - - -def put_payoff(s_t: float, strike: float, **kwargs): - r"""Computes the payoff of a European put option. - - Notes - ----- - - .. math:: - P(S_T, K) = \left(K-S_T, 0\right)^+ - - Parameters - ---------- - s_t : float - price - strike : float - the strike - - Returns - ------- - payoff : float - the payoff - """ - return np.maximum(strike - s_t, 0) - - -def futures_payoff(s_t: float, strike: float, **kwargs): - r"""Computes the payoff of a futures contract. - - Notes - ----- - - .. math:: - F(S_T, K) = \left(S_T-K, 0\right) - - Parameters - ---------- - s_t : float - price - strike : float - the strike - - Returns - ------- - payoff : float - the payoff - """ - return s_t - strike - - -def digital_call_payoff(s_t: float, strike: float, coupon: float = 1.0, **kwargs): - r"""Computes the payoff for a digital(binary) call option. - The formula is: - - Notes - ----- - - .. math:: - DC(S_T, K, Coupon) = Coupon\mathbb{1}_{\{S_T-K\}} - - Parameters - ---------- - s_t : float - current price of the underlying - strike : float - the strike - coupon : float - the amount received in case - that the underlying pring exceeds - the strike - - Returns - ------- - payoff : float - payoff of the european digital call option - """ - return np.where(s_t > strike, coupon, 0.0) - - -def digital_put_payoff(s_t: float, strike: float, coupon: float = 1.0, **kwargs): - r""" - Computes the payoff for a digital(binary) put option. - - Notes - ----- - The formula is: - - .. math:: - DC(S_T, K, Coupon) = Coupon\mathbb{1}_{\{K-S_T\}} - - Parameters - ---------- - s_t : float - current price of the underlying - strike : float - the strike - coupon : float - the amount received in case - that the underlying pring exceeds - the strike - - Returns - ------- - payoff : float - payoff of the european digital call option - """ - return np.where(s_t < strike, coupon, 0.0) - - -def bs_density( - s_t: float, - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - **kwargs -): - r""" - Evaluates the Black-Scholes density function at s_t for a given set - of parameters. The formula is: - - Notes - ----- - - .. math:: - \dfrac{1}{S_T\sigma\sqrt{2\pi T}}\exp\left(-\dfrac{ \ - \left(\log(S_T)-\mu\right)}{2\sigma^2T}\right) - - where :math:`\mu = (risk\_free\_rate-0.5\sigma)T+\log(S_0)`. - - Parameters - ---------- - s_t : float - point where we do the evaluation - s_0 : float - current price - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity: float - the maturity - - Returns - ------- - density : float - value of the Black-Scholes density function in s_t - """ - mean = (risk_free_rate - 0.5 * volatility * volatility) * maturity + np.log(s_0) - factor = s_t * volatility * np.sqrt(2 * np.pi * maturity) - exponent = -((np.log(s_t) - mean) ** 2) / (2 * volatility * volatility * maturity) - density = np.exp(exponent) / factor - return density - - -def bs_probability( - s_t: np.array, - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - **kwargs -): - r""" - Computes a discrete probability distribution from the Black-Scholes - density function for a given set of parameters. This is done by evaluating - the Black-Scholes density function in s_t and the normlising this result. - - Parameters - ---------- - s_t : numpy array - points where we define the discrete probability distribution - s_0 : float - current price - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity: float - the maturity - - Returns - ------- - distribution : numpy array - discrete probability distribution from Black-Scholes density - """ - density = bs_density(s_t, s_0, risk_free_rate, volatility, maturity) - return density / np.sum(density) - - -def bs_sde_solution( - x_: np.array, - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - **kwargs -): - r""" - For a certain parametrization :math:`x` it returns a value of the - underlying :math:`S_T(x)` and the probability density of that - value of the underlying. - - Notes - ----- - The formulas are: - - .. math:: - S_T = S_0e^{\sigma x+(risk_free_rate-\sigma^2/2)t} - p(S_T(x)) = N(x;mean = 0, variance = T) - - Parameters - ---------- - x_ : numpy array - parametrization - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity : float - the maturity - - Returns - ------- - s_t : numpy array - value of the underlying corresponding to parameter x - probability_density : numpy array - probability density corresponding to s_t - """ - probability = norm.pdf(x_) * np.sqrt(maturity) - probability = probability / np.sum(probability) - s_t = s_0 * np.exp( - volatility * x_ + (risk_free_rate - volatility * volatility / 2) * maturity - ) - return s_t, probability - - -def bs_call_price( - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - strike: float, - **kwargs -): - r""" - Computes the price for a European call option. - - Notes - ----- - - The formula is: - - .. math:: - C(S, T) = S\Phi(d_1)-Ke^{-rT}\Phi(d_2) - - Parameters - ---------- - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity : float - the maturity - strike : float - the strike - - Returns - ------- - price : float - price of the European call option - """ - first = np.log(s_0 / strike) - positive = (risk_free_rate + volatility * volatility / 2) * maturity - negative = (risk_free_rate - volatility * volatility / 2) * maturity - d_1 = (first + positive) / (volatility * np.sqrt(maturity)) - d_2 = (first + negative) / (volatility * np.sqrt(maturity)) - price = s_0 * norm.cdf(d_1) - strike * np.exp( - -risk_free_rate * maturity - ) * norm.cdf(d_2) - return price - - -def bs_put_price( - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - strike: float, - **kwargs -): - r""" - Computes the price for a European put option. - - Notes - ----- - The formula is: - - .. math:: - C(S, T) = Ke^{-rT}\Phi(-d_2)-S\Phi(-d_1) - - Parameters - ---------- - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity : float - the maturity - strike : float - the strike - - Returns - ------- - price : float - price of the European put option - """ - first = np.log(s_0 / strike) - positive = (risk_free_rate + volatility * volatility / 2) * maturity - negative = (risk_free_rate - volatility * volatility / 2) * maturity - d_1 = (first + positive) / (volatility * np.sqrt(maturity)) - d_2 = (first + negative) / (volatility * np.sqrt(maturity)) - price = strike * np.exp(-risk_free_rate * maturity) * norm.cdf( - -d_2 - ) - s_0 * norm.cdf(-d_1) - return price - - -def bs_digital_call_price( - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - strike: float, - coupon: float, - **kwargs -): - r""" - Computes the price for a digital(binary) call option. - - Notes - ----- - The formula is: - - .. math:: - DC(S, T) = e^{-rT}Coupon N(d_2) - - Parameters - ---------- - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity : float - the maturity - strike : float - the strike - coupon : float - the amount received in case that the underlying pring exceeds - the strike - - Returns - ------- - price : float - price of the digital call option - """ - first = np.log(s_0 / strike) - negative = (risk_free_rate - volatility * volatility / 2) * maturity - d_2 = (first + negative) / (volatility * np.sqrt(maturity)) - price = coupon * np.exp(-risk_free_rate * maturity) * norm.cdf(d_2) - return price - - -def bs_digital_put_price( - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - strike: float, - coupon: float, - **kwargs -): - r""" - Computes the price for a digital (binary) put option. - - Notes - ----- - The formula is: - - .. math:: - DC(S, T) = e^{-rT}Coupon N(-d_2) - - Parameters - ---------- - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity : float - the maturity - strike : float - the strike - coupon : float - the amount received in case that the underlying pring exceeds - the strike - - Returns - ------- - price : float - price of the digital call option - """ - first = np.log(s_0 / strike) - negative = (risk_free_rate - volatility * volatility / 2) * maturity - d_2 = (first + negative) / (volatility * np.sqrt(maturity)) - price = coupon * np.exp(-risk_free_rate * maturity) * norm.cdf(-d_2) - return price - - -def bs_exact_samples( - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - number_samples: int, - **kwargs -): - r""" - Computes samples from the exact solution of the Black-Scholes SDE. - - Notes - ----- - The formula is: - - .. math:: - S_T = S_0e^{\sigma*W_t+(risk_free_rate-\sigma^2/2)t} - - Parameters - ---------- - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity : float - the maturity - strike : float - the strike - number_samples : int - number of samples - - Returns - ------- - s_t : numpy array of floats - array of samples from the SDE. - """ - dw_t = np.random.randn(number_samples) * np.sqrt(maturity) - s_t = s_0 * np.exp( - volatility * dw_t + (risk_free_rate - volatility * volatility / 2) * maturity - ) - return s_t - - -def bs_em_samples( - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - number_samples: int, - time_steps: int, - **kwargs -): - r""" - Computes samples from the approximated solution of the Black-Scholes - SDE using the Euler-Maruyama discrimination. - - Notes - ----- - The formula is: - - .. math:: - S_{t+\Delta t} = S_t+rS_tdt+\sigma S_t N(0, 1)\sqrt{dt} - - Parameters - ---------- - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity : float - the maturity - strike : float - the strike - number_samples : int - number of samples - time steps : int - number of time steps - - Returns - ------- - s_t : numpy array of floats - array of samples from the SDE. - """ - dt = maturity / time_steps - s_t = np.ones(number_samples) * s_0 - for i in range(time_steps): - dw_t = np.random.randn(number_samples) * np.sqrt(dt) - s_t = s_t + risk_free_rate * s_t * dt + volatility * s_t * dw_t - return s_t - - -def bs_tree( - s_0: float, - risk_free_rate: float, - volatility: float, - maturity: float, - number_samples: int, - time_steps: int, - discretization: int, - bounds: float, - **kwargs -): - r""" - Computes the probabilities of all possible pahts from the - approximated solution of the Black-Scholes SDE using the - Euler-Maruyama discretization for a given discretization of the - Brownian motion. - - Parameters - ---------- - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - volatility : float - the volatility - maturity : float - the maturity - strike : float - the strike - number_samples : int - number of samples - time steps : int - number of time steps - discretization : float - number of points to build the discrete version - of the gaussian density - bounds : float - bounds of the Gaussian density - - Returns - ------- - s_t : numpy array of floats - array of samples from the SDE. - """ - dt = maturity / time_steps - x_ = np.linspace(-bounds, bounds, discretization) - p_x = norm.pdf(x_) - p_x = p_x / np.sum(p_x) - - s_t = [] - p_t = [] - s_t.append(np.array([s_0])) - p_t.append(np.array([1.0])) - for i in range(time_steps): - all_possible_paths = np.array(np.zeros(discretization ** (i + 1))) - all_possible_probabilities = np.array(np.zeros(discretization ** (i + 1))) - for j in range(len(s_t[i])): - single_possible_paths = ( - s_t[i][j] - + risk_free_rate * s_t[i][j] * dt - + volatility * s_t[i][j] * x_ * np.sqrt(dt) - ) - single_possible_probabilities = p_t[i][j] * p_x - - index = j * discretization - all_possible_paths[index : index + discretization] = single_possible_paths - all_possible_probabilities[ - index : index + discretization - ] = single_possible_probabilities - - s_t.append(all_possible_paths) - p_t.append(all_possible_probabilities) - return s_t, p_t - - -def geometric_sum(base: float, exponent: int, coeficient: float = 1.0, **kwargs): - r""" - Computes a geometric sum - - Notes - ----- - - .. math:: - s = a+a*base+a*base^2+a*base^3+...+a*base^{exponent} = \\ - = \sum_{k=0}^{exponent}a*base^k = \ - a\big(\frac{1-base^{exponent+1}}{1-base}\big) - - Parameters - ---------- - - base : float - base of the geometric sum - exponent: int - maximum exponent for the geometric sum - coeficient : float - coefficient for the geometric sum - - Returns - _______ - return : float - result of the geometric sum - - """ - return coeficient * (base ** (exponent + 1) - 1) / (base - 1) - - -def bs_forward_price( - s_0: float, risk_free_rate: float, maturity: float, strike: float, **kwargs -): - r""" - Computes the price for a forward contract. Note that it doesn't - assume that we are at the begining of the lifetime of the contract - - Notes - ----- - The formula is: - - .. math:: - V(S, T) = S-Ke^{-r(T-t)} - Parameters - ---------- - s_0 : float - current price of the underlying - risk_free_rate : float - risk free rate - maturity : float - time to expiration of the contract - strike : float - the strike - Returns - ------- - price : float - price of the forward - """ - price = s_0 - strike * np.exp(-risk_free_rate * maturity) - return price diff --git a/tnbs/BTC_02_AE/QQuantLib/finance/payoff_class.py b/tnbs/BTC_02_AE/QQuantLib/finance/payoff_class.py deleted file mode 100644 index 23d53c0..0000000 --- a/tnbs/BTC_02_AE/QQuantLib/finance/payoff_class.py +++ /dev/null @@ -1,101 +0,0 @@ -""" -Class definition for the PayOff - -Authors: Alberto Pedro Manzano Herrero & Gonzalo Ferro -""" -from functools import partial - -import QQuantLib.finance.classical_finance as cf -from QQuantLib.utils.utils import text_is_none - - -class PayOff: - - """ - Class for selecting derivative options and configuring them. - - Parameters - ---------- - - kwargs: dictionary - Dictionary for selecting and configuring the derivative option. \\ç - Implemented keys: - - * strike: float - strike of the option. - * coupon: float - only valid for Digital Options - * pay_off : string - type of pay_off function to load - """ - - def __init__(self, **kwargs): - """ - - Method for initializing the class - - """ - - self.pay_off_type = kwargs.get("pay_off_type", None) - text_is_none(self.pay_off_type, "pay_off_type", variable_type=str) - # European Options - self.strike = kwargs.get("strike", None) - # Digital Options - self.coupon = kwargs.get("coupon", None) - self.pay_off = None - self.pay_off_bs = None - self.get_pay_off() - - def get_pay_off(self): - """ - Select of a PayOff - """ - - if self.pay_off_type == "European_Call_Option": - text_is_none(self.strike, "strike", variable_type=float) - # PayOff - self.pay_off = partial(cf.call_payoff, strike=self.strike) - # Exact Solution for PayOff under BS - self.pay_off_bs = partial(cf.bs_call_price, strike=self.strike) - elif self.pay_off_type == "European_Put_Option": - text_is_none(self.strike, "strike", variable_type=float) - # PayOff - self.pay_off = partial(cf.put_payoff, strike=self.strike) - # Exact Solution for PayOff under BS - self.pay_off_bs = partial(cf.bs_put_price, strike=self.strike) - elif self.pay_off_type == "Digital_Call_Option": - text_is_none(self.strike, "strike", variable_type=float) - text_is_none(self.coupon, "coupon", variable_type=float) - # PayOff - self.pay_off = partial( - cf.digital_call_payoff, strike=self.strike, coupon=self.coupon - ) - # Exact Solution for PayOff under BS - self.pay_off_bs = partial( - cf.bs_digital_call_price, strike=self.strike, coupon=self.coupon - ) - elif self.pay_off_type == "Digital_Put_Option": - text_is_none(self.strike, "strike", variable_type=float) - text_is_none(self.coupon, "coupon", variable_type=float) - # PayOff - self.pay_off = partial( - cf.digital_put_payoff, strike=self.strike, coupon=self.coupon - ) - # Exact Solution for PayOff under BS - self.pay_off_bs = partial( - cf.bs_digital_put_price, strike=self.strike, coupon=self.coupon - ) - elif self.pay_off_type == "Futures": - text_is_none(self.strike, "strike", variable_type=float) - # PayOff - self.pay_off = partial(cf.futures_payoff, strike=self.strike) - # Exact Solution for PayOff under BS - self.pay_off_bs = partial(cf.bs_forward_price, strike=self.strike) - else: - raise ValueError( - "Not valid pay off type was provided. Valid types: \ - European_Call_Option,\n \ - European_Put_Option,\n \ - Digital_Call_Option,\n \ - Digital_Put_Option" - ) diff --git a/tnbs/BTC_02_AE/QQuantLib/finance/probability_class.py b/tnbs/BTC_02_AE/QQuantLib/finance/probability_class.py deleted file mode 100644 index bcdf333..0000000 --- a/tnbs/BTC_02_AE/QQuantLib/finance/probability_class.py +++ /dev/null @@ -1,144 +0,0 @@ -""" -Definition for DensityProbability Class. - -Authors: Alberto Pedro Manzano Herrero & Gonzalo Ferro -""" -from functools import partial - -import QQuantLib.finance.classical_finance as cf -from QQuantLib.utils.utils import text_is_none - - - -class DensityProbability: - - """ - Class for selecting pay off functions - algorithm - - Parameters - ---------- - - probability_type : string - type of probability density function to load - kwargs: dictionary - Dictionary for configuring the asset and the probability \\ - used for simulating its behaviour.Implemented keys: - - s_0 : float - initial value of the asset - risk_free_rate : float - risk free ratio - maturity : float - time where the probability wants to be calculated. - volatiliy : float - volatility of the asset. - """ - - def __init__(self, probability_type: str, **kwargs): - """ - - Method for initializing the class - - """ - - self.probability_type = probability_type - text_is_none(self.probability_type, "probability_type", variable_type=str) - self.probability = None - self.density_probability = None - self.probability = self.get_density( - self.probability_type, **kwargs) - self.density_probability = self.get_density_prob( - self.probability_type, **kwargs) - - @staticmethod - def get_density(probability_type, **kwargs): - """ - Create the probability function - - Parameters - ---------- - - probability_type : string - type of probability density function to load - kwargs: dictionary - with necessary information for configuring the probability \\ - density - - s_0 : float - initial value of the asset - risk_free_rate : float - risk free ratio - maturity : float - time where the probability wants to be calculated - volatiliy : float - volatility of the asset - - """ - - if probability_type == "Black-Scholes": - - s_0 = kwargs.get("s_0", None) - text_is_none(s_0, "s_0", variable_type=float) - risk_free_rate = kwargs.get("risk_free_rate", None) - text_is_none(risk_free_rate, "risk_free_rate", variable_type=float) - maturity = kwargs.get("maturity", None) - text_is_none(maturity, "maturity", variable_type=float) - volatility = kwargs.get("volatility", None) - text_is_none(volatility, "volatility", variable_type=float) - - else: - raise ValueError() - - return partial( - cf.bs_probability, - s_0=s_0, - risk_free_rate=risk_free_rate, - volatility=volatility, - maturity=maturity, - ) - - @staticmethod - def get_density_prob(probability_type, **kwargs): - """ - Configures a probability density - - Parameters - ---------- - - probability_type : string - type of probability density function to load - kwargs: dictionary - with necessary information for configuring the probability \\ - density. - - s_0 : float - initial value of the asset - risk_free_rate : float - risk free ratio - maturity : float - time where the probability wants to be calculated - volatiliy : float - volatility of the asset - """ - - if probability_type == "Black-Scholes": - - s_0 = kwargs.get("s_0", None) - text_is_none(s_0, "s_0", variable_type=float) - risk_free_rate = kwargs.get("risk_free_rate", None) - text_is_none(risk_free_rate, "risk_free_rate", variable_type=float) - maturity = kwargs.get("maturity", None) - text_is_none(maturity, "maturity", variable_type=float) - volatility = kwargs.get("volatility", None) - text_is_none(volatility, "volatility", variable_type=float) - - else: - raise ValueError() - return partial( - cf.bs_density, - s_0=s_0, - risk_free_rate=risk_free_rate, - volatility=volatility, - maturity=maturity, - ) diff --git a/tnbs/BTC_02_AE/QQuantLib/finance/quantum_integration.py b/tnbs/BTC_02_AE/QQuantLib/finance/quantum_integration.py index 663b3b1..84d5fb5 100644 --- a/tnbs/BTC_02_AE/QQuantLib/finance/quantum_integration.py +++ b/tnbs/BTC_02_AE/QQuantLib/finance/quantum_integration.py @@ -11,14 +11,10 @@ """ -import sys -import os -import re import warnings from copy import deepcopy import numpy as np import pandas as pd - from QQuantLib.DL.encoding_protocols import Encoding from QQuantLib.AE.ae_class import AE from QQuantLib.utils.utils import text_is_none @@ -49,7 +45,7 @@ def q_solve_integral(**kwargs): Note ---- - + Other kwargs input dictionary keys will be related with the encoding \\ of the integral into the quantum circuit \\ (see QQuantLib.DL.encoding_protocols) and for the configuration \\ @@ -67,9 +63,9 @@ def q_solve_integral(**kwargs): encoding = kwargs.get("encoding", None) ae_type = kwargs.get("ae_type", None) - if (encoding == 0) and (ae_type == "RQAE"): + if (encoding == 0) and (ae_type in ["RQAE", "eRQAE", "sRQAE", "mRQAE"]): string_error = ( - "RQAE method CAN NOT BE USED with encoding protocol: "+str(encoding) + "RQAEs methods CAN NOT BE USED with encoding protocol: "+str(encoding) ) warnings.warn(string_error) @@ -107,42 +103,57 @@ def q_solve_integral(**kwargs): linalg_qpu = kwargs.get("qpu", None) if linalg_qpu is None: raise ValueError("qpu is None. Please provide a valid qpu") - del kwargs['qpu'] + #del kwargs['qpu'] - ae_dict = deepcopy(kwargs) - ae_dict.update({"qpu": linalg_qpu}) + # ae_dict = deepcopy(kwargs) + #ae_dict.update({"qpu": linalg_qpu}) #Delete keys from encoding - for step in ["array_function", "array_probability", "encoding", "multiplexor"]: - ae_dict.pop(step, None) - ae_dict.pop("ae_type", None) + # for step in ["array_function", "array_probability", "encoding", "multiplexor"]: + # ae_dict.pop(step, None) + # ae_dict.pop("ae_type", None) #Instantiate AE solver solver_ae = AE( oracle=encode_class.oracle, target=encode_class.target, index=encode_class.index, - ae_type=ae_type, - **ae_dict) + #ae_type=ae_type, + **kwargs) # run the amplitude estimation algorithm solver_ae.run() # Recover amplitude estimation from ae_solver if encoding == 0: + # In the square encoding the np.sqrt(p_x * f_norm_x) is + # loaded in the amplitude of the state |x>^n|0>. + # The probability of measuring a |0> in the additional qubit + # is just sum_x(|p_x * f_norm_x|) ae_pdf = solver_ae.ae_pdf elif encoding == 1: - if ae_type == "RQAE": + if ae_type in ["RQAE", "eRQAE", "mRQAE", "sRQAE"]: #Amplitude is provided directly by this algorithm ae_pdf = solver_ae.ae_pdf else: #Other algorithms return probability ae_pdf = np.sqrt(solver_ae.ae_pdf) + if ae_type != "MLAE": + # We need to provided the mean of the square root bounds + # as the measured ae. Not valid for IQAE + ae_pdf["ae"] = (ae_pdf["ae_l"] + ae_pdf["ae_u"]) / 2.0 elif encoding == 2: - if ae_type == "RQAE": - #RQAE provides amplitude directly. + # In the direct encoding we load the sum_x(p_x * f_norm_x) + # into the amplitude of the |0>^{n+1} state + if ae_type in ["RQAE", "eRQAE", "mRQAE", "sRQAE"]: + # RQAE provides amplitude directly. ae_pdf = solver_ae.ae_pdf else: - #Other algorithms return probability + # AE algorithms provide an estimation of the probability + # so they provide: sum_x(p_x * f_norm_x) ** 2 + # we need to compute sqrt(sum_x(p_x * f_norm_x) ** 2) ae_pdf = np.sqrt(solver_ae.ae_pdf) + # We need to provided the mean of the square root bounds + # as the measured ae. + ae_pdf["ae"] = (ae_pdf["ae_l"] + ae_pdf["ae_u"]) / 2.0 else: raise ValueError("Not valid encoding key was provided!!!") #Now we need to deal with encoding normalisation diff --git a/tnbs/BTC_02_AE/QQuantLib/notebooks/01_AmplitudeEstimationKernel.ipynb b/tnbs/BTC_02_AE/QQuantLib/notebooks/01_AmplitudeEstimationKernel.ipynb index 2351818..a48cdcf 100644 --- a/tnbs/BTC_02_AE/QQuantLib/notebooks/01_AmplitudeEstimationKernel.ipynb +++ b/tnbs/BTC_02_AE/QQuantLib/notebooks/01_AmplitudeEstimationKernel.ipynb @@ -153,6 +153,14 @@ "* 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$**." ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f7bdbdbf", + "metadata": {}, + "outputs": [], + "source": [] } ], "metadata": { @@ -171,7 +179,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.9.9" } }, "nbformat": 4, diff --git a/tnbs/BTC_02_AE/QQuantLib/notebooks/02_AmplitudeEstimationBTC.ipynb b/tnbs/BTC_02_AE/QQuantLib/notebooks/02_AmplitudeEstimationBTC.ipynb index a7b278c..a7adb9e 100644 --- a/tnbs/BTC_02_AE/QQuantLib/notebooks/02_AmplitudeEstimationBTC.ipynb +++ b/tnbs/BTC_02_AE/QQuantLib/notebooks/02_AmplitudeEstimationBTC.ipynb @@ -496,7 +496,7 @@ "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 attributes of the encoding class:\n", "\n", "* *oracle*: myqlm gate with the unitary operator $\\mathbf{A(f_{x_i})}$\n", - "* *target*: the state $\\ket{\\Psi_0}$ of the unitary operator $\\mathbf{A(f_{x_i})}$\n", + "* *target*: the state $|\\Psi_0\\rangle$ of the unitary operator $\\mathbf{A(f_{x_i})}$\n", "* *index*: index affected by the unitary operator $\\mathbf{A(f_{x_i})}$\n", "\n", "So we need to provide the corresponding attributes of the object from the *Encoding* class to the **Grover** function for getting the desired **Grover** operator." @@ -546,7 +546,7 @@ "For each of the 2 $\\mathbf{A(f_{x_i})}$ operators and their correspondent Grover-like operators $\\mathbf{G}(\\mathbf{A(f_{x_i})})$ the desired *Amplitude Estimation* algorithm should be used. The Python class *AE* from *QQuantLib/AE/ae_class* module will be used for solving the **AE** problem. When this class is instantiated the following arguments should be provided:\n", "\n", "* *oracle*: *myqlm* gate with the unitary operator $\\mathbf{A(f_{x_i})}$\n", - "* *target*: the state $\\ket{\\Psi_0}$ of the unitary operator $\\mathbf{A(f_{x_i})}$\n", + "* *target*: the state $|\\Psi_0\\rangle$ of the unitary operator $\\mathbf{A(f_{x_i})}$\n", "* *index*: index affected by the unitary operator $\\mathbf{A(f_{x_i})}$\n", "* *ae_dictionary*: Python dictionary with the complete configuration of the **AE** algorithm. Most important key is: *ae_type* key where the **AE** algorithm should be specified. The following strings can be provided: **CQPEAE**, **IQPEAE**, **MLAE**, **IQAE**,**RQAE** and **MCAE**.\n", "\n", @@ -845,7 +845,8 @@ " 'ae_type': 'RQAE',\n", " 'epsilon' : 0.001,\n", " 'alpha': None,\n", - " 'gamma': 0.05\n", + " 'gamma': 0.05,\n", + " 'q': 2.0\n", "})\n", "ae_object = AE(\n", " oracle=encoding_object.oracle,\n", @@ -973,7 +974,8 @@ "ae_dictionary.update({\n", " 'ae_type': 'IQAE',\n", " 'epsilon' : 0.001,\n", - " 'alpha': 0.05\n", + " 'alpha': 0.05,\n", + " 'shots':100\n", "})\n", "\n", "encoding_dict = {\n", @@ -1207,7 +1209,8 @@ "ae_dictionary.update({\n", " 'ae_type': 'IQAE',\n", " 'epsilon' : 0.001,\n", - " 'alpha': 0.05\n", + " 'alpha': 0.05,\n", + " 'shots': 100\n", "})" ] }, @@ -1338,27 +1341,27 @@ "\n", "Examples: \n", "\n", - "For executing a 4 qubits discretization domain, computing the first interval integration using **MCAE** algorithm and the *c* QPU*:\n", + "For executing a 4 qubits discretization domain, computing the first interval integration using **MCAE** algorithm (with the configuration defined in ./jsons/integral_mcae_configuration.json) and the *c* QPU*:\n", "\n", - " python ae_sine_integral.py -n_qbits 4 -interval 0 -repetitions 1 -ae_type MCAE -qpu --exe\n", + " python ae_sine_integral.py -n_qbits 4 -interval 0 -repetitions 1 -qpu --exe -json_ae ./jsons/integral_mcae_configuration.json -id 0 --exe\n", "\n", - "For executing a 6 qubits discretization domain, computing the second interval integration using **RQAE** algorithm and the *c* QPU*:\n", + "For executing a 6 qubits discretization domain, computing the second interval integration using **RQAE** algorithm (with the configuration defined in ./jsons/integral_rqae_configuration.json) and the *c* QPU*:\n", "\n", - " python ae_sine_integral.py -n_qbits 6 -interval 1 -repetitions 1 -ae_type RQAE -qpu c --exe\n", + " python ae_sine_integral.py -n_qbits 6 -interval 1 -repetitions 1 -json_ae ./jsons/integral_rqae_configuration.json -qpu c --exe\n", "\n", "For configuring the **AE** algorithms you can edit the different JSON files under the folder **BTC_02_AE/jsons**. For each algorithm exists a JSON file:\n", "* integral_cqpeae_configuration.json: configure **CQPEAE**\n", "* integral_iqpeae_configuration.json: configure **IQPEAE**\n", - "* integral_mlae_configuration.json: configure **MLAE**\r", - "* integral_iqae_configuration.jso: configure **IQAE**\n", - "* integral_rqae_configuration.jso: configure **RQAE**\n", + "* integral_mlae_configuration.json: configure **MLAE**\n", + "* integral_iqae_configuration.json: configure **IQAE**\n", + "* integral_rqae_configuration.json: configure **RQAE**\n", "* integral_mcae_configuration.json: configure **MCAE**" ] }, { "cell_type": "code", "execution_count": null, - "id": "ac9d6c1c-fb7a-497f-9524-074bd0032879", + "id": "cd979893", "metadata": {}, "outputs": [], "source": [] @@ -1366,9 +1369,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:tnbs2c] *", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "conda-env-tnbs2c-py" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -1380,7 +1383,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.9.9" } }, "nbformat": 4, diff --git a/tnbs/BTC_02_AE/QQuantLib/notebooks/03_AmplitudeEstimationBenchmarkExecution.ipynb b/tnbs/BTC_02_AE/QQuantLib/notebooks/03_AmplitudeEstimationBenchmarkExecution.ipynb index 0b91192..bcfe56e 100644 --- a/tnbs/BTC_02_AE/QQuantLib/notebooks/03_AmplitudeEstimationBenchmarkExecution.ipynb +++ b/tnbs/BTC_02_AE/QQuantLib/notebooks/03_AmplitudeEstimationBenchmarkExecution.ipynb @@ -70,9 +70,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:tnbs] *", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "conda-env-tnbs-py" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -84,7 +84,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.9.9" } }, "nbformat": 4, diff --git a/tnbs/BTC_02_AE/QQuantLib/notebooks/04_AmpliutdeEstimationAlgorithms.ipynb b/tnbs/BTC_02_AE/QQuantLib/notebooks/04_AmpliutdeEstimationAlgorithms.ipynb index 088ff97..a3d3224 100644 --- a/tnbs/BTC_02_AE/QQuantLib/notebooks/04_AmpliutdeEstimationAlgorithms.ipynb +++ b/tnbs/BTC_02_AE/QQuantLib/notebooks/04_AmpliutdeEstimationAlgorithms.ipynb @@ -13,7 +13,7 @@ "\n", "$$|\\Psi\\rangle= \\mathbf{A}|0\\rangle_n = \\sqrt{a} |\\Psi_0\\rangle + \\sqrt{1-a}|\\Psi_1\\rangle \\tag{1}$$\n", "\n", - "The **AE kernel** tries to get an estimation of the probability of obtaining the state $\\ket{\\Psi_0}$, this is an estimator of $a$.\n", + "The **AE kernel** tries to get an estimation of the probability of obtaining the state $|\\Psi_0\\rangle$, this is an estimator of $a$.\n", "\n", "In the following cells we create the operator $\\mathbf{A(f_{x_i})}$ that we are going to use for testing the different algorithms." ] @@ -103,7 +103,7 @@ "source": [ "## 1. MonteCarlo Amplitude Estimation (MCAE)\n", "\n", - "In this case only the unitary operator $\\mathbf{A(f_{x_i})}$ that encodes the function to integrate is mandatory. the idea is to execute the circuit and measure the probability of obtaining the state $\\ket{\\Psi_0}$: $a$.\n", + "In this case only the unitary operator $\\mathbf{A(f_{x_i})}$ that encodes the function to integrate is mandatory. the idea is to execute the circuit and measure the probability of obtaining the state $|\\Psi_0\\rangle$: $a$.\n", "\n", "For this technique, the only mandatory input of the configuration dictionary of the **AE** class will be the number of shots for executing the circuit that is provided using the key: *shots*." ] @@ -243,7 +243,7 @@ "id": "f3b62422-2b09-4c60-aea1-4410907ff0ba", "metadata": {}, "source": [ - "## 3. Iterative Quantum Phase Estimation (IQAE).\n", + "## 3. Iterative Quantum Phase Estimation (IQPEAE).\n", "\n", "The **IQAE** is a variant of the **CQPEAE** algorithm where the Quantum Fourier Transform is done using only an additional qubit. A lower number of qubits are used but the depth of the circuits is much higher. This algorithm needs the Grover operator of $\\mathbf{A(f_{x_i})}$, $\\mathbf{G}(\\mathbf{A(f_{x_i})})$. The configuration keys for the *AE* are:\n", "\n", @@ -324,11 +324,11 @@ "\n", "This algorithm needs the Grover operator of $\\mathbf{A(f_{x_i})}$, $\\mathbf{G}(\\mathbf{A(f_{x_i})})$. It is an iterative algorithm. \n", "\n", - "In each step of the algorithm, the following operation is executed as a circuit: $$\\mathbf{G}^{m_k} \\mathbf{A(f_{x_i})} \\ket{0}_n$$ where $m_k$ is selected in advance and the circuit is measured a fixed number of shots ($n_k$) and the probability of obtaining the $\\ket{\\Psi_0}$ is obtained (h_k). \n", + "In each step of the algorithm, the following operation is executed as a circuit: $$\\mathbf{G}^{m_k} \\mathbf{A(f_{x_i})} |0\\rangle_n$$ where $m_k$ is selected in advance and the circuit is measured a fixed number of shots ($n_k$) and the probability of obtaining the $|\\Psi_0\\rangle$ is obtained (h_k). \n", "\n", "Using the properties of the *Grover* operator it can be shown that the following equation holds:\n", "\n", - "$$\\mathbf{G}^{m_k} \\mathbf{A(f_{x_i})} \\ket{0}_n =\\sin\\left((2m_k+1)\\theta\\right)|\\Psi_0\\rangle +\\cos\\left((2m_k+1)\\theta\\right)|\\Psi_1\\rangle,$$\n", + "$$\\mathbf{G}^{m_k} \\mathbf{A(f_{x_i})} |0\\rangle_n =\\sin\\left((2m_k+1)\\theta\\right)|\\Psi_0\\rangle +\\cos\\left((2m_k+1)\\theta\\right)|\\Psi_1\\rangle,$$\n", "\n", "where $$\\sqrt{a} = \\sin \\theta$$.\n", "\n", @@ -404,11 +404,11 @@ "id": "963d30c6-a89a-4664-b3a6-c08082cb0935", "metadata": {}, "source": [ - "## 5. Iterative Quantum Amplitude Estimation.\n", + "## 5. Iterative Quantum Amplitude Estimation (IQAE)\n", "\n", "This algorithm needs the Grover operator of $\\mathbf{A(f_{x_i})}$, $\\mathbf{G}(\\mathbf{A(f_{x_i})})$. It is an iterative algorithm.\n", "\n", - "In each step of the algorithm the property of the *Grover* operator: $$\\mathbf{G}^{m_k} \\mathbf{A(f_{x_i})} \\ket{0}_n =\\sin\\left((2m_k+1)\\theta\\right)|\\Psi_0\\rangle +\\cos\\left((2m_k+1)\\theta\\right)|\\Psi_1\\rangle,$$ is used for obtaining some close bounds for the estimation of $a$.\n", + "In each step of the algorithm the property of the *Grover* operator: $$\\mathbf{G}^{m_k} \\mathbf{A(f_{x_i})} |0\\rangle_n =\\sin\\left((2m_k+1)\\theta\\right)|\\Psi_0\\rangle +\\cos\\left((2m_k+1)\\theta\\right)|\\Psi_1\\rangle,$$ is used for obtaining some close bounds for the estimation of $a$.\n", "\n", "The **IQAE** needs an error $\\epsilon$ and a confident interval $\\alpha$ as input and the final mission is estimating some upper and lower bounds $(a_l, a_u)$ such $a$ satisfies that: $$P\\big[a \\in [a_l, a_u]\\big] \\gt 1-\\alpha$$ and $$\\frac{a_u-a_l}{2} \\leq \\epsilon$$\n", "\n", @@ -540,9 +540,9 @@ ], "metadata": { "kernelspec": { - "display_name": "Python [conda env:tnbs2c] *", + "display_name": "Python 3 (ipykernel)", "language": "python", - "name": "conda-env-tnbs2c-py" + "name": "python3" }, "language_info": { "codemirror_mode": { @@ -554,7 +554,7 @@ "name": "python", "nbconvert_exporter": "python", "pygments_lexer": "ipython3", - "version": "3.11.7" + "version": "3.9.9" } }, "nbformat": 4, diff --git a/tnbs/BTC_02_AE/QQuantLib/utils/data_extracting.py b/tnbs/BTC_02_AE/QQuantLib/utils/data_extracting.py index 5b2f651..6d6ffe7 100644 --- a/tnbs/BTC_02_AE/QQuantLib/utils/data_extracting.py +++ b/tnbs/BTC_02_AE/QQuantLib/utils/data_extracting.py @@ -6,15 +6,12 @@ Authors: Alberto Pedro Manzano Herrero & Gonzalo Ferro Costas """ -import sys -import os import time -import re +from copy import deepcopy import numpy as np import pandas as pd import qat.lang.AQASM as qlm from qat.core import Result - from QQuantLib.utils.utils import check_list_type pd.options.display.float_format = "{:.6f}".format @@ -58,7 +55,7 @@ def get_results( #print("BE AWARE!! linalg_qpu : {}".format(linalg_qpu)) # if type(quantum_object) == qlm.Program: if isinstance(quantum_object, qlm.Program): - q_prog = quantum_object + q_prog = deepcopy(quantum_object) arity = q_prog.qbit_count else: q_prog = create_qprogram(quantum_object) @@ -143,6 +140,7 @@ def create_qcircuit(prog_q): circuit : QLM circuit """ + #q_prog = deepcopy(prog_q) circuit = prog_q.to_circ(submatrices_only=True)#, inline=True) return circuit diff --git a/tnbs/BTC_02_AE/QQuantLib/utils/qlm_solver.py b/tnbs/BTC_02_AE/QQuantLib/utils/qlm_solver.py deleted file mode 100644 index 702e46d..0000000 --- a/tnbs/BTC_02_AE/QQuantLib/utils/qlm_solver.py +++ /dev/null @@ -1,57 +0,0 @@ -""" -This project has received funding from the European Union’s Horizon 2020 -research and innovation programme under Grant Agreement No. 951821 -https://www.neasqc.eu/ - -This module contains functions for calling QLM solver - -Authors: Alberto Pedro Manzano Herrero & Gonzalo Ferro Costas -""" - -from qat.qpus import get_default_qpu - - -def get_qpu(qpu=None): - """ - Function for selecting solver. - - Parameters - ---------- - - qpu : str - * qlmass: for trying to use QLM as a Service connection to CESGA QLM - * python: for using PyLinalg simulator. - * c: for using CLinalg simulator - - Returns - ---------- - - linal_qpu : solver for quantum jobs - """ - - if qpu is None: - raise ValueError( - "qpu CAN NOT BE NONE. Please select one of the three" + - " following options: qlmass, python, c") - if qpu == "qlmass": - try: - from qlmaas.qpus import LinAlg - linalg_qpu = LinAlg() - except (ImportError, OSError) as exception: - raise ImportError( - "Problem Using QLMaaS. Please create config file" + - "or use mylm solver") from exception - elif qpu == "python": - from qat.qpus import PyLinalg - linalg_qpu = PyLinalg() - elif qpu == "c": - from qat.qpus import CLinalg - linalg_qpu = CLinalg() - elif qpu == "default": - linalg_qpu = get_default_qpu() - else: - raise ValueError( - "Invalid value for qpu. Please select one of the three "+ - "following options: qlmass, python, c") - print("Following qpu will be used: {}".format(linalg_qpu)) - return linalg_qpu diff --git a/tnbs/BTC_02_AE/QQuantLib/utils/utils.py b/tnbs/BTC_02_AE/QQuantLib/utils/utils.py index 8e9142f..7c09ac6 100644 --- a/tnbs/BTC_02_AE/QQuantLib/utils/utils.py +++ b/tnbs/BTC_02_AE/QQuantLib/utils/utils.py @@ -2,10 +2,12 @@ This module contains several auxiliary functions needed by other scripts of the library + Authors: Alberto Pedro Manzano Herrero & Gonzalo Ferro Costas + Fast Walsh-Hadamard Transform is based on mex function written by Chengbo Li@Rice Uni for his TVAL3 algorithm: -https://github.com/dingluo/fwht/blob/master/FWHT.py + https://github.com/dingluo/fwht/blob/master/FWHT.py """ import time diff --git a/tnbs/BTC_02_AE/ae_sine_integral.py b/tnbs/BTC_02_AE/ae_sine_integral.py index e992f7b..ccabacf 100644 --- a/tnbs/BTC_02_AE/ae_sine_integral.py +++ b/tnbs/BTC_02_AE/ae_sine_integral.py @@ -38,9 +38,6 @@ def sine_integral(n_qbits, interval, ae_dictionary): DataFrame with the complete information of the benchmark """ - #Local copy for AE configuration dictionary - ae_dictionary_ = deepcopy(ae_dictionary) - start_time = time.time() #Section 2.1: Function for integration @@ -81,7 +78,7 @@ def sine_integral(n_qbits, interval, ae_dictionary): } #Now added the AE configuration. #The ae_dictionary_ has a local copy of the AE configuration. - q_solve_configuration.update(ae_dictionary_) + q_solve_configuration.update(ae_dictionary) #The q_solve_integral needs a QPU object. #q_solve_configuration["qpu"] = get_qpu(q_solve_configuration["qpu"]) @@ -107,7 +104,7 @@ def sine_integral(n_qbits, interval, ae_dictionary): #properly the KERNEL_BENCHMARK class #Adding the complete AE configuration - pdf = pd.DataFrame([ae_dictionary_]) + pdf = pd.DataFrame([ae_dictionary]) #Adding information about the computed integral pdf["interval"] = interval @@ -150,6 +147,70 @@ def sine_integral(n_qbits, interval, ae_dictionary): #elapsed_time, run_time, quantum_time] return pdf +def save(save, save_name, input_pdf, save_mode): + """ + For saving panda DataFrames to csvs + + Parameters + ---------- + + save: bool + For saving or not + save_nam: str + name for file + input_pdf: pandas DataFrame + save_mode: str + saving mode: overwrite (w) or append (a) + """ + if save: + with open(save_name, save_mode) as f_pointer: + input_pdf.to_csv( + f_pointer, + mode=save_mode, + header=f_pointer.tell() == 0, + sep=';' + ) + +def run_id( + n_qbits=None, + interval=0, + repetitions=None, + ae_config=None, + qpu=None, + save_=False, + folder_path=None, + ): + + ae_config.update({"qpu":get_qpu(qpu)}) + + if save_: + if folder_path is None: + raise ValueError("folder_name is None!") + if not os.path.exists(folder_path): + os.mkdir(folder_path) + + ae_type = ae_config["ae_type"] + base_name = ae_type + "_n_qbits_" + str(n_qbits) + \ + "_interval_" + str(interval) + ".csv" + file_name = folder_path + "/" + base_name + + list_of_pdfs = [] + for i in range(repetitions): + step_pdf = sine_integral( + n_qbits, + interval, + ae_config + ) + list_of_pdfs.append(step_pdf) + pdf = pd.concat(list_of_pdfs) + pdf.reset_index(drop=True, inplace=True) + save(save_, file_name, pdf, "w") + return pdf + + + + + def select_ae(ae_method): """ Function for selecting the AE algorithm used in the benchmark @@ -234,12 +295,18 @@ def select_ae(ae_method): ) #AE algorithm configuration arguments parser.add_argument( - "-ae_type", - dest="ae_type", + "-json_ae", + dest="json_ae", type=str, default=None, - help="AE algorithm for integral kernel: "+ - "[MLAE, IQAE, RQAE, MCAE, CQPEAE, IQPEAE]", + help="JSON AE algorithm configuration", + ) + parser.add_argument( + "--count", + dest="count", + default=False, + action="store_true", + help="For counting elements on the list", ) #For information about the configuation parser.add_argument( @@ -249,12 +316,19 @@ def select_ae(ae_method): action="store_true", help="For printing the AE algorihtm configuration." ) + parser.add_argument( + "-id", + dest="id", + type=int, + help="For executing only one element of the list", + default=None, + ) #QPU argument parser.add_argument( "-qpu", dest="qpu", type=str, - default="python", + default="c", help="QPU for simulation: See function get_qpu in get_qpu module", ) #Saving results arguments @@ -282,36 +356,55 @@ def select_ae(ae_method): ) args = parser.parse_args() - ae_configuration = select_ae(args.ae_type) - ae_configuration.update({"qpu":get_qpu(args.qpu)}) + with open(args.json_ae) as json_file: + ae_cfg = json.load(json_file) + # Creates the complete configuration for AE solvers + final_list = combination_for_list(ae_cfg) + if args.count: + print(len(final_list)) if args.print: - print(ae_configuration) + if args.id is not None: + print(final_list[args.id]) + else: + print(final_list) if args.execution: - list_of_pdfs = [] - for i in range(args.repetitions): - step_pdf = sine_integral( - args.n_qbits, - args.interval, - ae_configuration + if args.id is not None: + run_id( + n_qbits=args.n_qbits, + interval=args.interval, + repetitions=args.repetitions, + ae_config=final_list[args.id], + qpu=args.qpu, + folder_path=args.folder_path, + #qpu=args.qpu, + save_=args.save, ) - list_of_pdfs.append(step_pdf) - pdf = pd.concat(list_of_pdfs) - pdf.reset_index(drop=True, inplace=True) - print(pdf) - if args.save: - if args.folder_path is None: - raise ValueError("folder_name is None!") - if not os.path.exists(args.folder_path): - os.mkdir(args.folder_path) - base_name = args.ae_type + "_n_qbits_" + str(args.n_qbits) + \ - "_interval_" + str(args.interval) + ".csv" - file_name = args.folder_path + "/" + base_name - print(file_name) - with open(file_name, "w") as f_pointer: - pdf.to_csv( - f_pointer, - mode="w", - sep=';' - ) + # if args.execution: + # list_of_pdfs = [] + # for i in range(args.repetitions): + # step_pdf = sine_integral( + # args.n_qbits, + # args.interval, + # ae_configuration + # ) + # list_of_pdfs.append(step_pdf) + # pdf = pd.concat(list_of_pdfs) + # pdf.reset_index(drop=True, inplace=True) + # print(pdf) + # if args.save: + # if args.folder_path is None: + # raise ValueError("folder_name is None!") + # if not os.path.exists(args.folder_path): + # os.mkdir(args.folder_path) + # base_name = args.ae_type + "_n_qbits_" + str(args.n_qbits) + \ + # "_interval_" + str(args.interval) + ".csv" + # file_name = args.folder_path + "/" + base_name + # print(file_name) + # with open(file_name, "w") as f_pointer: + # pdf.to_csv( + # f_pointer, + # mode="w", + # sep=';' + # )