From be490aa1076a4c4fdbac9bc089910909b30caec6 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matthias=20H=C3=BCls?= Date: Tue, 5 Dec 2023 15:24:22 +0100 Subject: [PATCH 1/4] "test commit" --- src/quocslib/pulses/basis/PolynomialBasis.py | 89 ++++++++++++++++++++ src/quocslib/utils/map_dictionary.json | 3 + 2 files changed, 92 insertions(+) create mode 100644 src/quocslib/pulses/basis/PolynomialBasis.py diff --git a/src/quocslib/pulses/basis/PolynomialBasis.py b/src/quocslib/pulses/basis/PolynomialBasis.py new file mode 100644 index 0000000..fa762e5 --- /dev/null +++ b/src/quocslib/pulses/basis/PolynomialBasis.py @@ -0,0 +1,89 @@ +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Copyright 2021- QuOCS Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +import numpy as np + +from quocslib.pulses.BasePulse import BasePulse +from quocslib.pulses.basis.ChoppedBasis import ChoppedBasis +from quocslib.tools.randomgenerator import RandomNumberGenerator + + +class PolynomialBasis(ChoppedBasis): + """ + Dummy Basis class. It is used as a template for the creation of new basis. + """ + amplitude_variation: float + optimized_control_parameters: np.ndarray + optimized_super_parameters: np.ndarray + time_grid: np.ndarray + + def __init__(self, map_index: int, pulse_dictionary: dict, rng: RandomNumberGenerator = None, is_AD: bool = False): + """ + Constructor of the Polynomial Basis class. It calls the constructor of the parent class ChoppedBasis. + + :param int map_index: Index number to use to get the control parameter. + :param dict pulse_dictionary: The dictionary of the pulse is defined here. + :param RandomNumberGenerator rng: Random number generator. + :param bool is_AD: Flag to indicate if the pulse is used for the automatic differentiation. + """ + ################# + # Basis dependent settings + ################# + basis_dict = pulse_dictionary["basis"] + # Super Parameter number i.e. the basis vector number in the pulse parametrization + self.super_parameter_number = basis_dict.setdefault("basis_vector_number", 1) + # Number of control parameters to be optimized + self.control_parameters_number = 1 * self.super_parameter_number + ################# + # Standard Basis Settings: amplitude limits, amplitude variation for the simplex, + # distribution of super parameters, etc ... + ################ + # Constructor of the parent classes, i.e. Base Pulse and Chopped Basis + super().__init__(map_index=map_index, rng=rng, is_AD=is_AD, **pulse_dictionary) + ################# + # Basis dependent settings + ################# + # Scale coefficients: average distance of the points in the intial simplex + self.scale_coefficients = (self.amplitude_variation / np.sqrt(2) * np.ones((self.control_parameters_number, ))) + # Initial value of the parameters in the pulse parametrization + self.offset_coefficients = np.zeros((self.control_parameters_number, )) + + def _get_shaped_pulse(self) -> np.array: + """ + Definition of the pulse parametrization. It is called at every function evaluation to build the pulse and + return it as an array. + + :return np.array: The pulse as an array. + """ + ################# + # Standard Basis Settings: amplitude limits, amplitude variation for the simplex, + # distribution of super parameters, etc ... + ################ + # Pulse initialization + pulse = np.zeros(self.bins_number) + # Final time definition + final_time = self.final_time + # Pulse creation + xx = self.optimized_control_parameters + w = self.super_parameter_distribution_obj.w + t = self.time_grid + ################# + # Basis dependent settings + ################# + for ii in range(self.super_parameter_number): + pulse += xx[ii]*(t/final_time)**w[ii] + + return pulse \ No newline at end of file diff --git a/src/quocslib/utils/map_dictionary.json b/src/quocslib/utils/map_dictionary.json index a52d364..0c75816 100644 --- a/src/quocslib/utils/map_dictionary.json +++ b/src/quocslib/utils/map_dictionary.json @@ -34,6 +34,9 @@ "Sigmoid": {"module_name": "quocslib.pulses.basis.Sigmoid", "class_name": "Sigmoid"}, + "PolynomialBasis": + {"module_name": "quocslib.pulses.basis.PolynomialBasis", + "class_name": "PolynomialBasis"}, "Walsh": {"module_name": "quocslib.pulses.basis.Walsh", "class_name": "Walsh"} From 17a43e4b3f62900bd24c084ebc801be1b3dde897 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matthias=20H=C3=BCls?= Date: Thu, 7 Dec 2023 18:08:50 +0100 Subject: [PATCH 2/4] "just tests" --- .../Settings_in_Optimization_Dict.md | 2 +- .../IsingModelProblem.py | 185 ++++++++++++++++++ .../Some_Folder/Muh.py | 72 +++++++ .../Example_with_User_Basis/execute_dCRAB.py | 154 +++++++++++++++ .../settings_dCRAB.json | 61 ++++++ Examples/execute_dCRAB.py | 1 + Examples/testimporting.py | 9 + src/quocslib/utils/testimportfunction.py | 5 + 8 files changed, 488 insertions(+), 1 deletion(-) create mode 100644 Examples/Example_with_User_Basis/IsingModelProblem.py create mode 100644 Examples/Example_with_User_Basis/Some_Folder/Muh.py create mode 100644 Examples/Example_with_User_Basis/execute_dCRAB.py create mode 100644 Examples/Example_with_User_Basis/settings_dCRAB.json create mode 100644 Examples/testimporting.py create mode 100644 src/quocslib/utils/testimportfunction.py diff --git a/Documentation/Settings_in_Optimization_Dict.md b/Documentation/Settings_in_Optimization_Dict.md index 93d17a2..654d34d 100644 --- a/Documentation/Settings_in_Optimization_Dict.md +++ b/Documentation/Settings_in_Optimization_Dict.md @@ -15,7 +15,7 @@ Assuming you define the settings in the form of a .json file, the general struct "dump_format": "npz", # format of the results file "algorithm_settings": {...}, # settings related to the algorithm "pulses": [{...}, {...}, ...], # list of pulses and their settings - "parameters": [{...}, {...}, ...], # list of parameters and their settings + "parameters": [{...}, {...}, ...], # list of parameters and their settingsRaspberry Pi "times": [{...}, {...}, ...] # list of times and their settings } ~~~ diff --git a/Examples/Example_with_User_Basis/IsingModelProblem.py b/Examples/Example_with_User_Basis/IsingModelProblem.py new file mode 100644 index 0000000..34e2654 --- /dev/null +++ b/Examples/Example_with_User_Basis/IsingModelProblem.py @@ -0,0 +1,185 @@ +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Copyright 2021- QuOCS Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +import numpy as np +from quocslib.utils.AbstractFoM import AbstractFoM +from quocslib.timeevolution.piecewise_integrator import pw_evolution +from quocslib.tools.randomgenerator import RandomNumberGenerator +import functools + + +class IsingModel(AbstractFoM): + """A figure of merit class for optimization of the problem defined by Alastair Marshall via + https://arxiv.org/abs/2110.06187""" + + def __init__(self, args_dict: dict = None): + if args_dict is None: + args_dict = {} + ################################################################################################################ + # Dynamics variables + ################################################################################################################ + self.n_qubits = args_dict.setdefault("n_qubits", 5) + self.J = args_dict.setdefault("J", 1) + self.g = args_dict.setdefault("g", 2) + self.n_slices = args_dict.setdefault("n_slices", 100) + + self.H_drift = get_static_hamiltonian(self.n_qubits, self.J, self.g) + self.H_control = get_control_hamiltonian(self.n_qubits) + self.rho_0 = get_initial_state(self.n_qubits) + self.rho_target = get_target_state(self.n_qubits) + self.rho_final = np.zeros_like(self.rho_target) + # allocate a storage array + self.prop_store = [ + np.zeros_like(self.H_drift) for _ in range(self.n_slices) + ] + self.FoM_list = [] + self.rng = 0 + self.g_seed = args_dict.setdefault("g_seed", 0) + if self.g_seed != 0: + self.rng = RandomNumberGenerator(seed_number=self.g_seed) + self.g_variation = args_dict.setdefault("g_variation", 0) + self.stdev = args_dict.setdefault("stdev", 0.1) + + def get_control_Hamiltonians(self): + return self.H_control + + def get_drift_Hamiltonian(self): + if self.rng != 0: + return get_static_hamiltonian(self.n_qubits, self.J, + self.g + self.g_variation * (0.5 - self.rng.get_random_numbers(1)[0])) + else: + return get_static_hamiltonian(self.n_qubits, self.J, self.g) + + def get_target_state(self): + return self.rho_target + + def get_initial_state(self): + return self.rho_0 + + def get_propagator(self, + pulses_list: list = [], + time_grids_list: list = [], + parameters_list: list = []) -> np.array: + + drive = pulses_list[0].reshape(1, len(pulses_list[0])) + n_slices = self.n_slices + time_grid = time_grids_list[0] + # dt = time_grid[1] - time_grid[0] + dt = time_grid[-1] / len(time_grid) + # Compute the time evolution + self.prop_store = pw_evolution(self.prop_store, drive, + self.get_drift_Hamiltonian(), + [self.H_control], n_slices, dt) + return self.prop_store + + def get_FoM(self, + pulses: list = [], + parameters: list = [], + timegrids: list = []) -> dict: + """ """ + # Compute the final propagator + prop_store = self.get_propagator(pulses, timegrids, parameters) + U_final = functools.reduce(lambda a, b: a @ b, self.prop_store) + # evolve initial state + rho_final = U_final @ self.rho_0 @ U_final.T.conj() + # Calculate the fidelity + fidelity = fidelity_funct(rho_final.T, self.rho_target) + self.FoM_list.append(fidelity) + return {"FoM": -fidelity, "std": self.stdev} + + +i2 = np.eye(2) +sz = 0.5 * np.matrix([[1, 0], [0, -1]], dtype=np.complex128) +sx = 0.5 * np.matrix([[0, 1], [1, 0]], dtype=np.complex128) +psi0 = np.matrix([[1, 0], [0, 0]], dtype=np.complex128) +psiT = np.matrix([[0, 0], [0, 1]], dtype=np.complex128) + + +def tensor_together(A): + res = np.kron(A[0], A[1]) + if len(A) > 2: + for two in A[2:]: + res = np.kron(res, two) + else: + res = res + return res + + +def fidelity_funct(rho_evolved, rho_aim): + return np.abs(np.trace(rho_evolved.conj() @ rho_aim)) + + +def get_static_hamiltonian(nqu, J, g): + + dim = 2**nqu + H0 = np.zeros((dim, dim), dtype=np.complex128) + for j in range(nqu): + # set up holding array + rest = [i2] * nqu + # set the correct elements to sz + # check, so we can implement a loop around + if j == nqu - 1: + idx1 = j + idx2 = 0 + else: + idx1 = j + idx2 = j + 1 + rest[idx1] = sz + rest[idx2] = sz + H0 = H0 - J * tensor_together(rest) + + for j in range(nqu): + # set up holding array + rest = [i2] * nqu + # set the correct elements to sz + # check, so we can implement a loop around + if j == nqu - 1: + idx1 = j + idx2 = 1 + elif j == nqu - 2: + idx1 = j + idx2 = 0 + else: + idx1 = j + idx2 = j + 2 + rest[idx1] = sz + rest[idx2] = sz + H0 = H0 - g * tensor_together(rest) + return H0 + + +def get_control_hamiltonian(nqu: int): + # get the controls + dim = 2**nqu + H_at_t = np.zeros((dim, dim), dtype=np.complex128) + for j in range(nqu): + # set up holding array + rest = [i2] * nqu + # set the correct elements to sz + # check, so we can implement a loop around + rest[j] = sx + H_at_t = H_at_t + tensor_together(rest) + return H_at_t + + +def get_initial_state(nqu: int): + state = [psi0] * nqu + return tensor_together(state) + + +def get_target_state(nqu: int): + state = [psiT] * nqu + return tensor_together(state) diff --git a/Examples/Example_with_User_Basis/Some_Folder/Muh.py b/Examples/Example_with_User_Basis/Some_Folder/Muh.py new file mode 100644 index 0000000..32271a3 --- /dev/null +++ b/Examples/Example_with_User_Basis/Some_Folder/Muh.py @@ -0,0 +1,72 @@ +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Copyright 2021- QuOCS Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +import numpy as np + +from quocslib.pulses.basis.ChoppedBasis import ChoppedBasis +from quocslib.tools.randomgenerator import RandomNumberGenerator + + +class Muh(ChoppedBasis): + """ + Class for the Fourier basis. It inherits from the ChoppedBasis class. + """ + amplitude_variation: float + optimized_control_parameters: np.ndarray + optimized_super_parameters: np.ndarray + time_grid: np.ndarray + + def __init__(self, map_index: int, pulse_dictionary: dict, rng: RandomNumberGenerator = None, is_AD: bool = False): + """ + Constructor of the Fourier basis class. It calls the constructor of the parent class ChoppedBasis. + + :param int map_index: Index number to use to get the control parameter. + :param dict pulse_dictionary: The dictionary of the pulse is defined here. + :param RandomNumberGenerator rng: The random number generator. + :param bool is_AD: Flag to activate the automatic differentiation. + """ + basis_dict = pulse_dictionary["basis"] + # Frequencies number i.e. the basis vector number in the pulse parametrization + self.super_parameter_number = basis_dict.setdefault("basis_vector_number", 1) + # Number of control parameters to be optimized + self.control_parameters_number = 2 * self.super_parameter_number + # Constructor of the parent class, i.e. Chopped Basis + super().__init__(map_index=map_index, rng=rng, is_AD=is_AD, **pulse_dictionary) + # Define scale and offset coefficients + self.scale_coefficients = (self.amplitude_variation / np.sqrt(2) * np.ones((self.control_parameters_number,))) + self.offset_coefficients = np.zeros((self.control_parameters_number,)) + + print("HOOOWEEEYYY... THE Muh BASIS IS ACTUALLY BEING CALLED") + + def _get_shaped_pulse(self) -> np.array: + """ + Definition of the pulse parametrization. It is called at every function evaluation to build the pulse and + return it as an array. + + :return np.array: The pulse as an array. + """ + # Pulse definition + pulse = np.zeros(self.bins_number) + # Final time definition + final_time = self.final_time + # Pulse creation + xx = self.optimized_control_parameters + w = self.super_parameter_distribution_obj.w + t = self.time_grid + for ii in range(self.super_parameter_number): + pulse += xx[2 * ii] * np.sin(2 * np.pi * w[ii] * t / final_time) + xx[2 * ii + 1] * np.cos( + 2 * np.pi * w[ii] * t / final_time) + return pulse diff --git a/Examples/Example_with_User_Basis/execute_dCRAB.py b/Examples/Example_with_User_Basis/execute_dCRAB.py new file mode 100644 index 0000000..10facd8 --- /dev/null +++ b/Examples/Example_with_User_Basis/execute_dCRAB.py @@ -0,0 +1,154 @@ +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Copyright 2021- QuOCS Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +import os, sys, platform +import matplotlib.pyplot as plt +from quocslib.utils.inputoutput import readjson +from quocslib.Optimizer import Optimizer +from IsingModelProblem import IsingModel +import numpy as np +import time +import statistics + + +def plot_FoM(result_path, FoM_filename): + + if 'Windows' in platform.platform(): + opt_name = result_path.split('\\')[-1] + else: + opt_name = result_path.split('/')[-1] + + file_path = os.path.join(result_path, FoM_filename) + save_name = "FoM_" + opt_name + + FoM = [line.rstrip('\n') for line in open(file_path)] + FoM = [float(f) for f in FoM] + iterations = range(1, len(FoM) + 1) + # print('\nInitial FoM: %.4f' % FoM[0]) + # print('Final FoM: %.4f \n' % FoM[-1]) + min_FoM = min(FoM) + max_FoM = max(FoM) + difference = abs(max_FoM - min_FoM) + + fig = plt.figure(figsize=(11, 7)) + ax = fig.add_subplot(111) + plt.subplots_adjust(bottom=0.15, top=0.9, right=0.98, left=0.1) + + plt.plot(iterations, FoM, color='darkblue', linewidth=1.5, zorder=10) + # plt.scatter(x, y, color='k', s=15) + + plt.grid(True, which="both") + plt.ylim(min_FoM - 0.05 * difference, max_FoM + 0.05 * difference) + plt.xlabel('Iteration', fontsize=20) + plt.ylabel('FoM', fontsize=20) + # plt.savefig(os.path.join(folder, save_name + '.pdf')) + plt.savefig(os.path.join(result_path, save_name + '.png')) + + +def plot_controls(result_path): + + if 'Windows' in platform.platform(): + opt_name = result_path.split('\\')[-1] + else: + opt_name = result_path.split('/')[-1] + + for file in os.listdir(result_path): + if file.endswith('best_controls.npz'): + file_path = os.path.join(result_path, file) + + save_name = "Controls_" + opt_name + + controls = np.load(file_path) + + time_grid = [] + pulse = [] + + for data_name in controls.files: + if "time" in data_name: + time_grid = controls[data_name] + elif "pulse" in data_name: + pulse = controls[data_name] + + fig = plt.figure(figsize=(11, 7)) + ax = fig.add_subplot(111) + plt.subplots_adjust(bottom=0.15, top=0.9, right=0.98, left=0.1) + + plt.plot(time_grid, pulse, color='darkgreen', linewidth=1.5, zorder=10) + plt.grid(True, which="both") + plt.xlabel('Time', fontsize=20) + plt.ylabel('Amplitude', fontsize=20) + # plt.savefig(os.path.join(folder, save_name + '.pdf')) + plt.savefig(os.path.join(result_path, save_name + '.png')) + + +def main(optimization_dictionary: dict): + + args_dict = {"n_qubits": 5, "J": 1, "g": 2, "N_slices": 100, "T": 1.0, + "g_seed": 0, "g_variation": 0.1, "stdev": 0.01} + + optimization_dictionary["pulses"][0]["bins_number"] = args_dict["N_slices"] + optimization_dictionary["times"][0]["initial_value"] = args_dict["T"] + + if args_dict["g_seed"] != 0: + optimization_dictionary["algorithm_settings"]["re_evaluation"] = "{}" + + + # Create FoM object + FoM_object = IsingModel(args_dict=args_dict) + + # Define Optimizer + optimization_obj = Optimizer(optimization_dictionary, FoM_object) + + t1 = time.time() + + optimization_obj.execute() + + t2 = time.time() + + optimization_time = t2 - t1 + + with open( + os.path.join(optimization_obj.results_path, "optimization_time.txt"), "w" + ) as f: + f.write("# Time for optimization in seconds:\n") + f.write(str(optimization_time)) + + # fomlist = [element * (-1) for element in optimization_obj.fom_list] + fomlist = [element for element in FoM_object.FoM_list] + np.savetxt(os.path.join(optimization_obj.results_path, "fom.txt"), fomlist) + + plot_FoM(optimization_obj.results_path, "fom.txt") + # plot_controls(optimization_obj.results_path) + + opt_controls = optimization_obj.opt_alg_obj.get_best_controls() + + statistics_fom_list = [] + num_for_average = 50 + for i in range(num_for_average): + statistics_fom_list.append(FoM_object.get_FoM(**opt_controls)["FoM"]*(-1)) + + mittel = statistics.mean(statistics_fom_list) + deviation = statistics.stdev(statistics_fom_list) + + with open(os.path.join(optimization_obj.results_path, "statistics.txt"), 'w') as f: + f.write('averaged over {} evals\n'.format(num_for_average)) + f.write('mean:{}\n'.format(mittel)) + f.write('stdev: {}\n'.format(deviation)) + + +if __name__ == "__main__": + main(readjson(os.path.join(os.getcwd(), "settings_dCRAB.json"))) + diff --git a/Examples/Example_with_User_Basis/settings_dCRAB.json b/Examples/Example_with_User_Basis/settings_dCRAB.json new file mode 100644 index 0000000..da651dd --- /dev/null +++ b/Examples/Example_with_User_Basis/settings_dCRAB.json @@ -0,0 +1,61 @@ +{ + "optimization_client_name": "Optimization_dCRAB_IsingModel", + "create_logfile": false, + "algorithm_settings": { + "algorithm_name": "dCRAB", + "super_iteration_number": 3, + "max_eval_total": 100, + "dsm_settings": { + "general_settings": { + "dsm_algorithm_name": "NelderMead", + "is_adaptive": true + }, + "stopping_criteria": { + "xatol": 1e-14, + "frtol": 1e-3, + "change_based_stop": { + "cbs_funct_evals": 200, + "cbs_change": 0.01 + }, + "max_eval": 1000 + } + } + }, + "pulses": [ + { + "pulse_name": "Pulse1", + "upper_limit": 1000.0, + "lower_limit": -1000.0, + "time_name": "time1", + "amplitude_variation": 10.0, + "basis": { + "basis_module": "Some_Folder.Muh", + "basis_class": "Muh", + "basis_name": "Muh", + "basis_vector_number": 5, + "random_super_parameter_distribution": { + "distribution_name": "Uniform", + "lower_limit": 0.01, + "upper_limit": 5.0 + } + }, + "scaling_function": { + "function_type": "lambda_function", + "lambda_function": "lambda t: 1.0 + 0.0*t" + }, + "initial_guess": { + "function_type": "lambda_function", + "lambda_function": "lambda t: 0.0 + 0.0*t" + } + } + ], + "times": [ + { + "time_name": "time1" + } + ], + "parameters": [], + "communication": { + "communication_type": "AllInOneCommunication" + } +} \ No newline at end of file diff --git a/Examples/execute_dCRAB.py b/Examples/execute_dCRAB.py index 0cc2429..558163c 100644 --- a/Examples/execute_dCRAB.py +++ b/Examples/execute_dCRAB.py @@ -148,6 +148,7 @@ def main(optimization_dictionary: dict): print("\nBest FoM: {}".format(optimization_obj.opt_alg_obj.best_FoM)) +print(__name__) if __name__ == "__main__": main(readjson(os.path.join(os.getcwd(), "settings_dCRAB.json"))) diff --git a/Examples/testimporting.py b/Examples/testimporting.py new file mode 100644 index 0000000..4bf4c68 --- /dev/null +++ b/Examples/testimporting.py @@ -0,0 +1,9 @@ +from matplotlib import pyplot as plt +import numpy as np + +from quocslib.utils.testimportfunction import gaussian + +xs = np.linspace(0,1,100) + +plt.plot(xs, gaussian(xs,1,0.5,0)) +plt.show() \ No newline at end of file diff --git a/src/quocslib/utils/testimportfunction.py b/src/quocslib/utils/testimportfunction.py new file mode 100644 index 0000000..987cd6b --- /dev/null +++ b/src/quocslib/utils/testimportfunction.py @@ -0,0 +1,5 @@ +import numpy as np + + +def gaussian(x,a,b,c): + return a * np.exp(-x**2/(b**2)) + c \ No newline at end of file From a739728ac11a78de592eba34d799368967b0d01f Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matthias=20H=C3=BCls?= Date: Thu, 14 Dec 2023 14:32:49 +0100 Subject: [PATCH 3/4] Implemented a continuation option --- .../Settings_in_Optimization_Dict.md | 3 + Examples/execute_dCRAB.py | 2 +- Examples/testimporting.py | 9 - src/quocslib/Optimizer.py | 9 +- .../communication/AllInOneCommunication.py | 29 ++- src/quocslib/tools/logger.py | 3 +- src/quocslib/utils/Import_previous_results.py | 57 +++++ src/quocslib/utils/testimportfunction.py | 5 - tests/test_continuation.py | 231 ++++++++++++++++++ 9 files changed, 326 insertions(+), 22 deletions(-) delete mode 100644 Examples/testimporting.py create mode 100644 src/quocslib/utils/Import_previous_results.py delete mode 100644 src/quocslib/utils/testimportfunction.py create mode 100644 tests/test_continuation.py diff --git a/Documentation/Settings_in_Optimization_Dict.md b/Documentation/Settings_in_Optimization_Dict.md index 654d34d..601790e 100644 --- a/Documentation/Settings_in_Optimization_Dict.md +++ b/Documentation/Settings_in_Optimization_Dict.md @@ -13,6 +13,7 @@ Assuming you define the settings in the form of a .json file, the general struct "create_logfile": true, # determines if you want to save the log-file "console_info": true, # determines if you want the optimization output printed to the console "dump_format": "npz", # format of the results file + "continuation_datetime": "YYYYmmdd_HHMMSS" # date of optimization to be continued "algorithm_settings": {...}, # settings related to the algorithm "pulses": [{...}, {...}, ...], # list of pulses and their settings "parameters": [{...}, {...}, ...], # list of parameters and their settingsRaspberry Pi @@ -30,6 +31,8 @@ The `console_info` key determines if the optimization output is shown in the ter The `dump_format` key specifies the format of the results file (best controls and some meta data). Currently you can choose between "npz" and "json". The default (if you do not give this key) is "npz". +The `continuation_datetime` key determines wether an existing optimization job should be continued and the results saved in the corresponding folder. If intended, the corresponding date and time should be provided in the format "YYYYmmdd_HHMMSS" and the corresponding job name speficied in `optimization_client_name`. QuOCS will then try to import the best controls found in the previous run as initial values for the continuation. Here make sure, to use the same names for pulses and parameters in the `optimiaztion_dictionary` as in the previous run. The default of this key is "no" and will create a new folder with a timestamp of the point of creation. + **Tip:** You can also change specific entries in the code after reading in the .json file if you, e.g., want to sweep certain parameters or have the name of the optimization defined on runtime. diff --git a/Examples/execute_dCRAB.py b/Examples/execute_dCRAB.py index 558163c..08a903b 100644 --- a/Examples/execute_dCRAB.py +++ b/Examples/execute_dCRAB.py @@ -148,7 +148,7 @@ def main(optimization_dictionary: dict): print("\nBest FoM: {}".format(optimization_obj.opt_alg_obj.best_FoM)) -print(__name__) + if __name__ == "__main__": main(readjson(os.path.join(os.getcwd(), "settings_dCRAB.json"))) diff --git a/Examples/testimporting.py b/Examples/testimporting.py deleted file mode 100644 index 4bf4c68..0000000 --- a/Examples/testimporting.py +++ /dev/null @@ -1,9 +0,0 @@ -from matplotlib import pyplot as plt -import numpy as np - -from quocslib.utils.testimportfunction import gaussian - -xs = np.linspace(0,1,100) - -plt.plot(xs, gaussian(xs,1,0.5,0)) -plt.show() \ No newline at end of file diff --git a/src/quocslib/Optimizer.py b/src/quocslib/Optimizer.py index 3529fed..165073c 100644 --- a/src/quocslib/Optimizer.py +++ b/src/quocslib/Optimizer.py @@ -21,6 +21,7 @@ from quocslib.communication.AllInOneCommunication import AllInOneCommunication from quocslib.utils.BestDump import BestDump from quocslib.utils.AbstractFoM import AbstractFoM +from quocslib.utils.Import_previous_results import update_opti_dict class Optimizer: @@ -58,6 +59,8 @@ def __init__(self, self.dump_format = optimization_dict.setdefault("dump_format", "npz") self.optimization_direction = optimization_dict["algorithm_settings"].setdefault("optimization_direction", "minimization") + self.continuation_datetime = optimization_dict.setdefault("continuation_datetime", "no") + self.communication_obj = AllInOneCommunication(interface_job_name=self.interface_job_name, FoM_obj=FoM_object, handle_exit_obj=handle_exit_obj, @@ -66,7 +69,11 @@ def __init__(self, create_logfile=self.create_logfile, console_info=self.console_info, dump_format=self.dump_format, - optimization_direction=self.optimization_direction) + optimization_direction=self.optimization_direction, + continuation_datetime=self.continuation_datetime) + + if self.communication_obj.is_continuation: + optimization_dict = update_opti_dict(optimization_dict, self.communication_obj) self.results_path = self.communication_obj.results_path diff --git a/src/quocslib/communication/AllInOneCommunication.py b/src/quocslib/communication/AllInOneCommunication.py index 93e4d40..69e0b14 100644 --- a/src/quocslib/communication/AllInOneCommunication.py +++ b/src/quocslib/communication/AllInOneCommunication.py @@ -35,7 +35,8 @@ def __init__(self, create_logfile: bool = True, console_info: bool = True, dump_format: str = "npz", - optimization_direction: str = "minimization"): + optimization_direction: str = "minimization", + continuation_datetime: str = "no"): """ In case the user chooses to run the optimization in his device, this class is used by the OptimizationAlgorithm. The objects to dump the results, calculate the figure of merit, and the logger are created here. @@ -57,25 +58,43 @@ def __init__(self, (self.message_signal, self.FoM_plot_signal, self.controls_update_signal) = comm_signals_list # Pre job name pre_job_name = interface_job_name + # Optimization folder name + optimization_folder = "QuOCS_Results" # Datetime for 1-1 association - self.date_time = str(time.strftime("%Y%m%d_%H%M%S")) + # Check, if optimization is continuation + self.is_continuation = False + queued_logger_info = None + if continuation_datetime == "no": + self.date_time = str(time.strftime("%Y%m%d_%H%M%S")) + else: + continuation_job_folder = os.path.join(os.getcwd(), optimization_folder, continuation_datetime + "_" + pre_job_name) + if not os.path.isdir(continuation_job_folder): + queued_logger_info = "Continuation attempt: No Folder " + continuation_job_folder + " found, new folder created" + self.date_time = str(time.strftime("%Y%m%d_%H%M%S")) + else: + self.date_time = continuation_datetime + self.is_continuation = True + queued_logger_info = "Continue optimization from " + continuation_datetime + "_" + pre_job_name # Client job name to send to the Server self.client_job_name = self.date_time + "_" + pre_job_name ### # Logging, Results, Figure of merit evaluation ... ### - # Optimization folder - optimization_folder = "QuOCS_Results" + # Optimization folder self.results_path = os.path.join(os.getcwd(), optimization_folder, self.client_job_name) if not os.path.isdir(os.path.join(os.getcwd(), optimization_folder)): os.makedirs(os.path.join(os.getcwd(), optimization_folder)) # Create the folder for logging and results - os.makedirs(self.results_path) + if not os.path.isdir(self.results_path): + os.makedirs(self.results_path) # Write the current quocs lib version in the file with open(os.path.join(self.results_path, "quocs_version.txt"), "w") as version_file: version_file.write("QuOCS library version: {0}".format(quocslib_version)) # Create logging object self.logger = create_logger(self.results_path, self.date_time, create_logfile=create_logfile, console_info=console_info) + # print queued logger info + if not queued_logger_info == None: + self.logger.info(queued_logger_info) # Print function evaluation and figure of merit self.print_general_log = True # Figure of merit object diff --git a/src/quocslib/tools/logger.py b/src/quocslib/tools/logger.py index d39ff1b..0f6356f 100644 --- a/src/quocslib/tools/logger.py +++ b/src/quocslib/tools/logger.py @@ -52,7 +52,8 @@ def create_logger(results_path, date_time, create_logfile=True, console_info=Tru console_handler.setFormatter(logging.Formatter(print_format)) # Log file handler if create_logfile: - file_handler = logging.FileHandler(log_filename) + # file_handler = logging.FileHandler(log_filename) + file_handler = logging.FileHandler(log_filename, mode='a') file_handler.setLevel(logging.INFO) file_handler.setFormatter(logging.Formatter(log_format, date_format)) # Add handler for logfile to the logger diff --git a/src/quocslib/utils/Import_previous_results.py b/src/quocslib/utils/Import_previous_results.py new file mode 100644 index 0000000..91a2aa8 --- /dev/null +++ b/src/quocslib/utils/Import_previous_results.py @@ -0,0 +1,57 @@ +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Copyright 2021- QuOCS Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +import os +import numpy as np +from quocslib.utils.inputoutput import readjson + + +def update_opti_dict(optimization_dict: dict, comm_obj: object) -> dict: + """ + Load optimal results of previous run and initial guess into an optimization dictionary + :param dict optimization_dictionary: optimization_dictionary to be updated + :comm_obj: communication object of optimization + :return dict: updated optimization dictionary + """ + + try: + if optimization_dict["dump_format"] == "json": + best_res_path = os.path.join(comm_obj.results_path, comm_obj.date_time + "_best_controls.json") + best_res = readjson(best_res_path) + else: + best_res_path = os.path.join(comm_obj.results_path, comm_obj.date_time + "_best_controls.npz") + best_res = np.load(best_res_path) + + for pulse in optimization_dict["pulses"]: # use same optimization dictionary as in previous optimization + pulse_name = pulse["pulse_name"] # make sure to use the same pulse names in opti_dict as in best_controls + prev_opt_pulse = best_res[pulse_name] + initial_guess = {"function_type": "list_function", "list_function": prev_opt_pulse} + pulse["initial_guess"] = initial_guess + comm_obj.logger.info(f"Initial guess for pulse {pulse_name} imported from previous results") + + for param in optimization_dict["parameters"]: + param_name = param["parameter_name"] + prev_opt_param = best_res[param_name] + param["initial_value"] = prev_opt_param + comm_obj.logger.info(f"Initial guess for parameter {param_name} imported from previous results") + + except: + comm_obj.logger.warn("Previous optimal controls could not be imported for continuation." + " Check if ..._best_controls...-file exists and pulse/paremeter names coincide with optimization dictionary") + + return optimization_dict + + diff --git a/src/quocslib/utils/testimportfunction.py b/src/quocslib/utils/testimportfunction.py deleted file mode 100644 index 987cd6b..0000000 --- a/src/quocslib/utils/testimportfunction.py +++ /dev/null @@ -1,5 +0,0 @@ -import numpy as np - - -def gaussian(x,a,b,c): - return a * np.exp(-x**2/(b**2)) + c \ No newline at end of file diff --git a/tests/test_continuation.py b/tests/test_continuation.py new file mode 100644 index 0000000..2aae363 --- /dev/null +++ b/tests/test_continuation.py @@ -0,0 +1,231 @@ +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ +# Copyright 2021- QuOCS Team +# +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ + +import os, platform +import matplotlib.pyplot as plt +import numpy as np +from quocslib.optimalcontrolproblems.OneQubitProblem import OneQubit +from quocslib.optimalcontrolproblems.IsingModelProblem import IsingModel +from quocslib.Optimizer import Optimizer +import pytest + + +def plot_FoM(result_path, FoM): + + save_name = "FoM_plot" + + iterations = range(1, len(FoM) + 1) + min_FoM = min(FoM) + max_FoM = max(FoM) + difference = abs(max_FoM - min_FoM) + + fig = plt.figure(figsize=(11, 7)) + ax = fig.add_subplot(111) + plt.subplots_adjust(bottom=0.15, top=0.9, right=0.98, left=0.1) + + plt.plot(iterations, FoM, color='darkblue', linewidth=1.5, zorder=10) + + plt.grid(True, which="both") + plt.ylim(min_FoM - 0.05 * difference, max_FoM + 0.05 * difference) + plt.xlabel('Iteration', fontsize=20) + plt.ylabel('FoM', fontsize=20) + plt.savefig(os.path.join(result_path, save_name + '.png')) + +def run_dCRAB_opti(optimization_dictionary): + # define some parameters for the optimization + args_dict = { + "initial_state": "[1.0 , 0.0]", + "target_state": "[1.0/np.sqrt(2), -1j/np.sqrt(2)]", + "optimization_factor": -1.0 + } + # Create FoM object + FoM_object = OneQubit(args_dict=args_dict) + + # Define Optimizer + optimization_obj = Optimizer(optimization_dictionary, FoM_object) + + initial_guess = optimization_obj.get_optimization_algorithm().controls.pulse_objs_list[0].initial_guess_pulse + initial_param = optimization_obj.get_optimization_algorithm().controls.parameter_objs_list[0].value + + # Run optimization + + optimization_obj.execute() + + optimal_pulse = optimization_obj.get_optimization_algorithm().get_best_controls()["pulses"][0] + optimal_param = optimization_obj.get_optimization_algorithm().get_best_controls()["parameters"][0] + + fomlist = optimization_obj.get_optimization_algorithm().FoM_list + + res_path = optimization_obj.results_path + datetime = optimization_obj.communication_obj.date_time + + return res_path, datetime, initial_guess, optimal_pulse, initial_param, optimal_param, fomlist + + + +def test_dCRAB_continuation(): + + optimization_dictionary = { + "optimization_client_name": "continuation_test", + "optimization_direction": "minimization", + "continuation_datetime": "no", + "dump_format": "json", + "algorithm_settings": { + "algorithm_name": "dCRAB", + "super_iteration_number": 2, + "max_eval_total": 100, + "FoM_goal": 0.00001, + "dsm_settings": { + "general_settings": { + "dsm_algorithm_name": "NelderMead", + "is_adaptive": True + }, + "stopping_criteria": { + "max_eval": 50, + } + }, + "random_number_generator": { + "seed_number": 42 + } + }, + "pulses": [{ + "pulse_name": "Pulse_1", + "upper_limit": 5.0, + "lower_limit": -5.0, + "bins_number": 101, + "time_name": "time_1", + "amplitude_variation": 5.0, + "basis": { + "basis_name": "Fourier", + "basis_vector_number": 2, + "random_super_parameter_distribution": { + "distribution_name": "Uniform", + "lower_limit": 0.1, + "upper_limit": 5.0 + } + }, + "initial_guess": { + "function_type": "lambda_function", + "lambda_function": "lambda t: np.pi/3.0 + 0.0*t" + } + }], + "parameters": [{"parameter_name": "Parameter0", + "lower_limit": -2.0, + "upper_limit": 2.0, + "initial_value": 0.4, + "amplitude_variation": 0.5}], + "times": [{ + "time_name": "time_1", + "initial_value": 3.0 + }] + } + res_path1, datetime1, _ , optimal_pulse1, _, optimal_param1, fomlist1 = run_dCRAB_opti(optimization_dictionary) + + optimization_dictionary["continuation_datetime"] = datetime1 + + # optimization_dictionary["pulses"][0]["pulse_name"] = "new_name" + + res_path2, datetime2, inital_guess2, optimal_pulse2, initial_param2, optimal_param2, fomlist2 = run_dCRAB_opti(optimization_dictionary) + + res_path3, datetime3, inital_guess3 , _ ,initial_param3, _ , fomlist3 = run_dCRAB_opti(optimization_dictionary) + + plot_FoM(res_path2, fomlist1 + fomlist2 + fomlist3) + + + assert res_path1 == res_path2 == res_path3 # test, if similar result path for both optimizations + assert datetime1 == datetime2 == datetime3 # test for similar datetime + assert np.array_equal(inital_guess2, optimal_pulse1) # test, if optimal pulse is given imported as initial guess + assert np.array_equal(inital_guess3, optimal_pulse2) + assert initial_param2 == optimal_param1 + assert initial_param3 == optimal_param2 + assert min(fomlist3) <= min(fomlist2) <= min(fomlist1) # test, if second optimization improved the results + + + + +def run_GRAPE_opti(optimization_dictionary): + + FoM_object = IsingModel(args_dict={}) + + optimization_obj = Optimizer(optimization_dictionary, FoM_object) + + initial_guess = optimization_obj.get_optimization_algorithm().controls.pulse_objs_list[0].initial_guess_pulse + # Run optimization + + optimization_obj.execute() + + optimal_pulse = optimization_obj.get_optimization_algorithm().get_best_controls()["pulses"][0] + fomlist = optimization_obj.get_optimization_algorithm().FoM_list + + res_path = optimization_obj.results_path + datetime = optimization_obj.communication_obj.date_time + + return res_path, datetime, initial_guess, optimal_pulse, fomlist + + +def test_GRAPE_continuation(): + optimization_dictionary = { + "optimization_client_name": "continuation_test", + "optimization_direction": "minimization", + "continuation_datetime": "no", + "dump_format": "json", + "algorithm_settings": { + "algorithm_name": "GRAPE", + "stopping_criteria": {"max_eval_total": 150} + }, + "pulses": [{ + "pulse_name": "Pulse_1", + "upper_limit": 100.0, + "lower_limit": -100.0, + "bins_number": 100, + "amplitude_variation": 20.0, + "time_name": "time_1", + "basis": { + "basis_name": "PiecewiseBasis", + "bins_number": 100 + }, + "initial_guess": { + "function_type": "lambda_function", + "lambda_function": "lambda t: 0.0 + 0.0*t" + } + }], + "parameters": [], + "times": [{ + "time_name": "time_1", + "initial_value": 1.0 + }] + } + + res_path1, datetime1, _ , optimal_pulse1, fomlist1 = run_GRAPE_opti(optimization_dictionary) + + optimization_dictionary["continuation_datetime"] = datetime1 + + # optimization_dictionary["pulses"][0]["pulse_name"] = "new_name" + + res_path2, datetime2, inital_guess2, optimal_pulse2, fomlist2 = run_GRAPE_opti(optimization_dictionary) + + res_path3, datetime3, inital_guess3 , _ , fomlist3 = run_GRAPE_opti(optimization_dictionary) + + plot_FoM(res_path2, fomlist1 + fomlist2 + fomlist3) + + + assert res_path1 == res_path2 == res_path3 # test, if similar result path for both optimizations + assert datetime1 == datetime2 == datetime3 # test for similar datetime + assert np.array_equal(inital_guess2, optimal_pulse1) # test, if optimal pulse is given imported as initial guess + assert np.array_equal(inital_guess3, optimal_pulse2) + assert min(fomlist3) <= min(fomlist2) <= min(fomlist1) # test, if second optimization improved the results + +test_dCRAB_continuation() \ No newline at end of file From 73cdfc7e454b0347cca09b744c25e66fe1ae5146 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Matthias=20H=C3=BCls?= Date: Thu, 14 Dec 2023 15:26:26 +0100 Subject: [PATCH 4/4] corrected previous changes --- .../Settings_in_Optimization_Dict.md | 2 +- .../IsingModelProblem.py | 185 ------------------ .../Some_Folder/Muh.py | 72 ------- .../Example_with_User_Basis/execute_dCRAB.py | 154 --------------- .../settings_dCRAB.json | 61 ------ src/quocslib/pulses/basis/PolynomialBasis.py | 89 --------- src/quocslib/utils/map_dictionary.json | 3 - 7 files changed, 1 insertion(+), 565 deletions(-) delete mode 100644 Examples/Example_with_User_Basis/IsingModelProblem.py delete mode 100644 Examples/Example_with_User_Basis/Some_Folder/Muh.py delete mode 100644 Examples/Example_with_User_Basis/execute_dCRAB.py delete mode 100644 Examples/Example_with_User_Basis/settings_dCRAB.json delete mode 100644 src/quocslib/pulses/basis/PolynomialBasis.py diff --git a/Documentation/Settings_in_Optimization_Dict.md b/Documentation/Settings_in_Optimization_Dict.md index 601790e..484e25b 100644 --- a/Documentation/Settings_in_Optimization_Dict.md +++ b/Documentation/Settings_in_Optimization_Dict.md @@ -16,7 +16,7 @@ Assuming you define the settings in the form of a .json file, the general struct "continuation_datetime": "YYYYmmdd_HHMMSS" # date of optimization to be continued "algorithm_settings": {...}, # settings related to the algorithm "pulses": [{...}, {...}, ...], # list of pulses and their settings - "parameters": [{...}, {...}, ...], # list of parameters and their settingsRaspberry Pi + "parameters": [{...}, {...}, ...], # list of parameters and their settings "times": [{...}, {...}, ...] # list of times and their settings } ~~~ diff --git a/Examples/Example_with_User_Basis/IsingModelProblem.py b/Examples/Example_with_User_Basis/IsingModelProblem.py deleted file mode 100644 index 34e2654..0000000 --- a/Examples/Example_with_User_Basis/IsingModelProblem.py +++ /dev/null @@ -1,185 +0,0 @@ -# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -# Copyright 2021- QuOCS Team -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -import numpy as np -from quocslib.utils.AbstractFoM import AbstractFoM -from quocslib.timeevolution.piecewise_integrator import pw_evolution -from quocslib.tools.randomgenerator import RandomNumberGenerator -import functools - - -class IsingModel(AbstractFoM): - """A figure of merit class for optimization of the problem defined by Alastair Marshall via - https://arxiv.org/abs/2110.06187""" - - def __init__(self, args_dict: dict = None): - if args_dict is None: - args_dict = {} - ################################################################################################################ - # Dynamics variables - ################################################################################################################ - self.n_qubits = args_dict.setdefault("n_qubits", 5) - self.J = args_dict.setdefault("J", 1) - self.g = args_dict.setdefault("g", 2) - self.n_slices = args_dict.setdefault("n_slices", 100) - - self.H_drift = get_static_hamiltonian(self.n_qubits, self.J, self.g) - self.H_control = get_control_hamiltonian(self.n_qubits) - self.rho_0 = get_initial_state(self.n_qubits) - self.rho_target = get_target_state(self.n_qubits) - self.rho_final = np.zeros_like(self.rho_target) - # allocate a storage array - self.prop_store = [ - np.zeros_like(self.H_drift) for _ in range(self.n_slices) - ] - self.FoM_list = [] - self.rng = 0 - self.g_seed = args_dict.setdefault("g_seed", 0) - if self.g_seed != 0: - self.rng = RandomNumberGenerator(seed_number=self.g_seed) - self.g_variation = args_dict.setdefault("g_variation", 0) - self.stdev = args_dict.setdefault("stdev", 0.1) - - def get_control_Hamiltonians(self): - return self.H_control - - def get_drift_Hamiltonian(self): - if self.rng != 0: - return get_static_hamiltonian(self.n_qubits, self.J, - self.g + self.g_variation * (0.5 - self.rng.get_random_numbers(1)[0])) - else: - return get_static_hamiltonian(self.n_qubits, self.J, self.g) - - def get_target_state(self): - return self.rho_target - - def get_initial_state(self): - return self.rho_0 - - def get_propagator(self, - pulses_list: list = [], - time_grids_list: list = [], - parameters_list: list = []) -> np.array: - - drive = pulses_list[0].reshape(1, len(pulses_list[0])) - n_slices = self.n_slices - time_grid = time_grids_list[0] - # dt = time_grid[1] - time_grid[0] - dt = time_grid[-1] / len(time_grid) - # Compute the time evolution - self.prop_store = pw_evolution(self.prop_store, drive, - self.get_drift_Hamiltonian(), - [self.H_control], n_slices, dt) - return self.prop_store - - def get_FoM(self, - pulses: list = [], - parameters: list = [], - timegrids: list = []) -> dict: - """ """ - # Compute the final propagator - prop_store = self.get_propagator(pulses, timegrids, parameters) - U_final = functools.reduce(lambda a, b: a @ b, self.prop_store) - # evolve initial state - rho_final = U_final @ self.rho_0 @ U_final.T.conj() - # Calculate the fidelity - fidelity = fidelity_funct(rho_final.T, self.rho_target) - self.FoM_list.append(fidelity) - return {"FoM": -fidelity, "std": self.stdev} - - -i2 = np.eye(2) -sz = 0.5 * np.matrix([[1, 0], [0, -1]], dtype=np.complex128) -sx = 0.5 * np.matrix([[0, 1], [1, 0]], dtype=np.complex128) -psi0 = np.matrix([[1, 0], [0, 0]], dtype=np.complex128) -psiT = np.matrix([[0, 0], [0, 1]], dtype=np.complex128) - - -def tensor_together(A): - res = np.kron(A[0], A[1]) - if len(A) > 2: - for two in A[2:]: - res = np.kron(res, two) - else: - res = res - return res - - -def fidelity_funct(rho_evolved, rho_aim): - return np.abs(np.trace(rho_evolved.conj() @ rho_aim)) - - -def get_static_hamiltonian(nqu, J, g): - - dim = 2**nqu - H0 = np.zeros((dim, dim), dtype=np.complex128) - for j in range(nqu): - # set up holding array - rest = [i2] * nqu - # set the correct elements to sz - # check, so we can implement a loop around - if j == nqu - 1: - idx1 = j - idx2 = 0 - else: - idx1 = j - idx2 = j + 1 - rest[idx1] = sz - rest[idx2] = sz - H0 = H0 - J * tensor_together(rest) - - for j in range(nqu): - # set up holding array - rest = [i2] * nqu - # set the correct elements to sz - # check, so we can implement a loop around - if j == nqu - 1: - idx1 = j - idx2 = 1 - elif j == nqu - 2: - idx1 = j - idx2 = 0 - else: - idx1 = j - idx2 = j + 2 - rest[idx1] = sz - rest[idx2] = sz - H0 = H0 - g * tensor_together(rest) - return H0 - - -def get_control_hamiltonian(nqu: int): - # get the controls - dim = 2**nqu - H_at_t = np.zeros((dim, dim), dtype=np.complex128) - for j in range(nqu): - # set up holding array - rest = [i2] * nqu - # set the correct elements to sz - # check, so we can implement a loop around - rest[j] = sx - H_at_t = H_at_t + tensor_together(rest) - return H_at_t - - -def get_initial_state(nqu: int): - state = [psi0] * nqu - return tensor_together(state) - - -def get_target_state(nqu: int): - state = [psiT] * nqu - return tensor_together(state) diff --git a/Examples/Example_with_User_Basis/Some_Folder/Muh.py b/Examples/Example_with_User_Basis/Some_Folder/Muh.py deleted file mode 100644 index 32271a3..0000000 --- a/Examples/Example_with_User_Basis/Some_Folder/Muh.py +++ /dev/null @@ -1,72 +0,0 @@ -# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -# Copyright 2021- QuOCS Team -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -import numpy as np - -from quocslib.pulses.basis.ChoppedBasis import ChoppedBasis -from quocslib.tools.randomgenerator import RandomNumberGenerator - - -class Muh(ChoppedBasis): - """ - Class for the Fourier basis. It inherits from the ChoppedBasis class. - """ - amplitude_variation: float - optimized_control_parameters: np.ndarray - optimized_super_parameters: np.ndarray - time_grid: np.ndarray - - def __init__(self, map_index: int, pulse_dictionary: dict, rng: RandomNumberGenerator = None, is_AD: bool = False): - """ - Constructor of the Fourier basis class. It calls the constructor of the parent class ChoppedBasis. - - :param int map_index: Index number to use to get the control parameter. - :param dict pulse_dictionary: The dictionary of the pulse is defined here. - :param RandomNumberGenerator rng: The random number generator. - :param bool is_AD: Flag to activate the automatic differentiation. - """ - basis_dict = pulse_dictionary["basis"] - # Frequencies number i.e. the basis vector number in the pulse parametrization - self.super_parameter_number = basis_dict.setdefault("basis_vector_number", 1) - # Number of control parameters to be optimized - self.control_parameters_number = 2 * self.super_parameter_number - # Constructor of the parent class, i.e. Chopped Basis - super().__init__(map_index=map_index, rng=rng, is_AD=is_AD, **pulse_dictionary) - # Define scale and offset coefficients - self.scale_coefficients = (self.amplitude_variation / np.sqrt(2) * np.ones((self.control_parameters_number,))) - self.offset_coefficients = np.zeros((self.control_parameters_number,)) - - print("HOOOWEEEYYY... THE Muh BASIS IS ACTUALLY BEING CALLED") - - def _get_shaped_pulse(self) -> np.array: - """ - Definition of the pulse parametrization. It is called at every function evaluation to build the pulse and - return it as an array. - - :return np.array: The pulse as an array. - """ - # Pulse definition - pulse = np.zeros(self.bins_number) - # Final time definition - final_time = self.final_time - # Pulse creation - xx = self.optimized_control_parameters - w = self.super_parameter_distribution_obj.w - t = self.time_grid - for ii in range(self.super_parameter_number): - pulse += xx[2 * ii] * np.sin(2 * np.pi * w[ii] * t / final_time) + xx[2 * ii + 1] * np.cos( - 2 * np.pi * w[ii] * t / final_time) - return pulse diff --git a/Examples/Example_with_User_Basis/execute_dCRAB.py b/Examples/Example_with_User_Basis/execute_dCRAB.py deleted file mode 100644 index 10facd8..0000000 --- a/Examples/Example_with_User_Basis/execute_dCRAB.py +++ /dev/null @@ -1,154 +0,0 @@ -# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -# Copyright 2021- QuOCS Team -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -import os, sys, platform -import matplotlib.pyplot as plt -from quocslib.utils.inputoutput import readjson -from quocslib.Optimizer import Optimizer -from IsingModelProblem import IsingModel -import numpy as np -import time -import statistics - - -def plot_FoM(result_path, FoM_filename): - - if 'Windows' in platform.platform(): - opt_name = result_path.split('\\')[-1] - else: - opt_name = result_path.split('/')[-1] - - file_path = os.path.join(result_path, FoM_filename) - save_name = "FoM_" + opt_name - - FoM = [line.rstrip('\n') for line in open(file_path)] - FoM = [float(f) for f in FoM] - iterations = range(1, len(FoM) + 1) - # print('\nInitial FoM: %.4f' % FoM[0]) - # print('Final FoM: %.4f \n' % FoM[-1]) - min_FoM = min(FoM) - max_FoM = max(FoM) - difference = abs(max_FoM - min_FoM) - - fig = plt.figure(figsize=(11, 7)) - ax = fig.add_subplot(111) - plt.subplots_adjust(bottom=0.15, top=0.9, right=0.98, left=0.1) - - plt.plot(iterations, FoM, color='darkblue', linewidth=1.5, zorder=10) - # plt.scatter(x, y, color='k', s=15) - - plt.grid(True, which="both") - plt.ylim(min_FoM - 0.05 * difference, max_FoM + 0.05 * difference) - plt.xlabel('Iteration', fontsize=20) - plt.ylabel('FoM', fontsize=20) - # plt.savefig(os.path.join(folder, save_name + '.pdf')) - plt.savefig(os.path.join(result_path, save_name + '.png')) - - -def plot_controls(result_path): - - if 'Windows' in platform.platform(): - opt_name = result_path.split('\\')[-1] - else: - opt_name = result_path.split('/')[-1] - - for file in os.listdir(result_path): - if file.endswith('best_controls.npz'): - file_path = os.path.join(result_path, file) - - save_name = "Controls_" + opt_name - - controls = np.load(file_path) - - time_grid = [] - pulse = [] - - for data_name in controls.files: - if "time" in data_name: - time_grid = controls[data_name] - elif "pulse" in data_name: - pulse = controls[data_name] - - fig = plt.figure(figsize=(11, 7)) - ax = fig.add_subplot(111) - plt.subplots_adjust(bottom=0.15, top=0.9, right=0.98, left=0.1) - - plt.plot(time_grid, pulse, color='darkgreen', linewidth=1.5, zorder=10) - plt.grid(True, which="both") - plt.xlabel('Time', fontsize=20) - plt.ylabel('Amplitude', fontsize=20) - # plt.savefig(os.path.join(folder, save_name + '.pdf')) - plt.savefig(os.path.join(result_path, save_name + '.png')) - - -def main(optimization_dictionary: dict): - - args_dict = {"n_qubits": 5, "J": 1, "g": 2, "N_slices": 100, "T": 1.0, - "g_seed": 0, "g_variation": 0.1, "stdev": 0.01} - - optimization_dictionary["pulses"][0]["bins_number"] = args_dict["N_slices"] - optimization_dictionary["times"][0]["initial_value"] = args_dict["T"] - - if args_dict["g_seed"] != 0: - optimization_dictionary["algorithm_settings"]["re_evaluation"] = "{}" - - - # Create FoM object - FoM_object = IsingModel(args_dict=args_dict) - - # Define Optimizer - optimization_obj = Optimizer(optimization_dictionary, FoM_object) - - t1 = time.time() - - optimization_obj.execute() - - t2 = time.time() - - optimization_time = t2 - t1 - - with open( - os.path.join(optimization_obj.results_path, "optimization_time.txt"), "w" - ) as f: - f.write("# Time for optimization in seconds:\n") - f.write(str(optimization_time)) - - # fomlist = [element * (-1) for element in optimization_obj.fom_list] - fomlist = [element for element in FoM_object.FoM_list] - np.savetxt(os.path.join(optimization_obj.results_path, "fom.txt"), fomlist) - - plot_FoM(optimization_obj.results_path, "fom.txt") - # plot_controls(optimization_obj.results_path) - - opt_controls = optimization_obj.opt_alg_obj.get_best_controls() - - statistics_fom_list = [] - num_for_average = 50 - for i in range(num_for_average): - statistics_fom_list.append(FoM_object.get_FoM(**opt_controls)["FoM"]*(-1)) - - mittel = statistics.mean(statistics_fom_list) - deviation = statistics.stdev(statistics_fom_list) - - with open(os.path.join(optimization_obj.results_path, "statistics.txt"), 'w') as f: - f.write('averaged over {} evals\n'.format(num_for_average)) - f.write('mean:{}\n'.format(mittel)) - f.write('stdev: {}\n'.format(deviation)) - - -if __name__ == "__main__": - main(readjson(os.path.join(os.getcwd(), "settings_dCRAB.json"))) - diff --git a/Examples/Example_with_User_Basis/settings_dCRAB.json b/Examples/Example_with_User_Basis/settings_dCRAB.json deleted file mode 100644 index da651dd..0000000 --- a/Examples/Example_with_User_Basis/settings_dCRAB.json +++ /dev/null @@ -1,61 +0,0 @@ -{ - "optimization_client_name": "Optimization_dCRAB_IsingModel", - "create_logfile": false, - "algorithm_settings": { - "algorithm_name": "dCRAB", - "super_iteration_number": 3, - "max_eval_total": 100, - "dsm_settings": { - "general_settings": { - "dsm_algorithm_name": "NelderMead", - "is_adaptive": true - }, - "stopping_criteria": { - "xatol": 1e-14, - "frtol": 1e-3, - "change_based_stop": { - "cbs_funct_evals": 200, - "cbs_change": 0.01 - }, - "max_eval": 1000 - } - } - }, - "pulses": [ - { - "pulse_name": "Pulse1", - "upper_limit": 1000.0, - "lower_limit": -1000.0, - "time_name": "time1", - "amplitude_variation": 10.0, - "basis": { - "basis_module": "Some_Folder.Muh", - "basis_class": "Muh", - "basis_name": "Muh", - "basis_vector_number": 5, - "random_super_parameter_distribution": { - "distribution_name": "Uniform", - "lower_limit": 0.01, - "upper_limit": 5.0 - } - }, - "scaling_function": { - "function_type": "lambda_function", - "lambda_function": "lambda t: 1.0 + 0.0*t" - }, - "initial_guess": { - "function_type": "lambda_function", - "lambda_function": "lambda t: 0.0 + 0.0*t" - } - } - ], - "times": [ - { - "time_name": "time1" - } - ], - "parameters": [], - "communication": { - "communication_type": "AllInOneCommunication" - } -} \ No newline at end of file diff --git a/src/quocslib/pulses/basis/PolynomialBasis.py b/src/quocslib/pulses/basis/PolynomialBasis.py deleted file mode 100644 index fa762e5..0000000 --- a/src/quocslib/pulses/basis/PolynomialBasis.py +++ /dev/null @@ -1,89 +0,0 @@ -# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ -# Copyright 2021- QuOCS Team -# -# Licensed under the Apache License, Version 2.0 (the "License"); -# you may not use this file except in compliance with the License. -# You may obtain a copy of the License at -# -# http://www.apache.org/licenses/LICENSE-2.0 -# -# Unless required by applicable law or agreed to in writing, software -# distributed under the License is distributed on an "AS IS" BASIS, -# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. -# See the License for the specific language governing permissions and -# limitations under the License. -# ++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++++ - -import numpy as np - -from quocslib.pulses.BasePulse import BasePulse -from quocslib.pulses.basis.ChoppedBasis import ChoppedBasis -from quocslib.tools.randomgenerator import RandomNumberGenerator - - -class PolynomialBasis(ChoppedBasis): - """ - Dummy Basis class. It is used as a template for the creation of new basis. - """ - amplitude_variation: float - optimized_control_parameters: np.ndarray - optimized_super_parameters: np.ndarray - time_grid: np.ndarray - - def __init__(self, map_index: int, pulse_dictionary: dict, rng: RandomNumberGenerator = None, is_AD: bool = False): - """ - Constructor of the Polynomial Basis class. It calls the constructor of the parent class ChoppedBasis. - - :param int map_index: Index number to use to get the control parameter. - :param dict pulse_dictionary: The dictionary of the pulse is defined here. - :param RandomNumberGenerator rng: Random number generator. - :param bool is_AD: Flag to indicate if the pulse is used for the automatic differentiation. - """ - ################# - # Basis dependent settings - ################# - basis_dict = pulse_dictionary["basis"] - # Super Parameter number i.e. the basis vector number in the pulse parametrization - self.super_parameter_number = basis_dict.setdefault("basis_vector_number", 1) - # Number of control parameters to be optimized - self.control_parameters_number = 1 * self.super_parameter_number - ################# - # Standard Basis Settings: amplitude limits, amplitude variation for the simplex, - # distribution of super parameters, etc ... - ################ - # Constructor of the parent classes, i.e. Base Pulse and Chopped Basis - super().__init__(map_index=map_index, rng=rng, is_AD=is_AD, **pulse_dictionary) - ################# - # Basis dependent settings - ################# - # Scale coefficients: average distance of the points in the intial simplex - self.scale_coefficients = (self.amplitude_variation / np.sqrt(2) * np.ones((self.control_parameters_number, ))) - # Initial value of the parameters in the pulse parametrization - self.offset_coefficients = np.zeros((self.control_parameters_number, )) - - def _get_shaped_pulse(self) -> np.array: - """ - Definition of the pulse parametrization. It is called at every function evaluation to build the pulse and - return it as an array. - - :return np.array: The pulse as an array. - """ - ################# - # Standard Basis Settings: amplitude limits, amplitude variation for the simplex, - # distribution of super parameters, etc ... - ################ - # Pulse initialization - pulse = np.zeros(self.bins_number) - # Final time definition - final_time = self.final_time - # Pulse creation - xx = self.optimized_control_parameters - w = self.super_parameter_distribution_obj.w - t = self.time_grid - ################# - # Basis dependent settings - ################# - for ii in range(self.super_parameter_number): - pulse += xx[ii]*(t/final_time)**w[ii] - - return pulse \ No newline at end of file diff --git a/src/quocslib/utils/map_dictionary.json b/src/quocslib/utils/map_dictionary.json index 0c75816..a52d364 100644 --- a/src/quocslib/utils/map_dictionary.json +++ b/src/quocslib/utils/map_dictionary.json @@ -34,9 +34,6 @@ "Sigmoid": {"module_name": "quocslib.pulses.basis.Sigmoid", "class_name": "Sigmoid"}, - "PolynomialBasis": - {"module_name": "quocslib.pulses.basis.PolynomialBasis", - "class_name": "PolynomialBasis"}, "Walsh": {"module_name": "quocslib.pulses.basis.Walsh", "class_name": "Walsh"}