diff --git a/tests/interactive/CRAB_gate_closed.md b/tests/interactive/CRAB_gate_closed.md index ea36556..221b42b 100644 --- a/tests/interactive/CRAB_gate_closed.md +++ b/tests/interactive/CRAB_gate_closed.md @@ -5,14 +5,14 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.17.1 + jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# CRAB algorithm for a closed system +# CRAB algorithm for a closed system (gate synthesis) ```python import matplotlib.pyplot as plt @@ -94,7 +94,7 @@ plt.show() ```python assert res_crab.infidelity < 0.001 -assert fidelity(evolution.states[-1], target_gate) > 1-0.001 +assert np.allclose(fidelity(evolution.states[-1], target_gate), 1 - res_crab.infidelity, atol=1e-3) ``` ```python diff --git a/tests/interactive/CRAB_gate_open.md b/tests/interactive/CRAB_gate_open.md new file mode 100644 index 0000000..b31f27c --- /dev/null +++ b/tests/interactive/CRAB_gate_open.md @@ -0,0 +1,115 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.17.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# CRAB algorithm for an open system (gate synthesis) + +```python +import matplotlib.pyplot as plt +import numpy as np +from qutip import gates, qeye, liouvillian, sigmam, sigmax, sigmay, sigmaz +import qutip as qt +from qutip_qoc import Objective, optimize_pulses + +def fidelity(gate_super, target_super): + gate_oper = qt.Qobj(gate_super.data) + target_oper = qt.Qobj(target_super.data) + + return np.abs(gate_oper.overlap(target_oper) / target_oper.norm()) +``` + +## Problem setup + + +```python +omega = 0.1 # energy splitting +gamma = 0.1 # amplitude damping +sx, sy, sz = sigmax(), sigmay(), sigmaz() +c_ops = [np.sqrt(gamma) * sigmam()] + +Hd = 1 / 2 * omega * sz +Hd = liouvillian(H=Hd, c_ops=c_ops) +Hc = [liouvillian(sx), liouvillian(sy), liouvillian(sz)] +H = [Hd, Hc[0], Hc[1], Hc[2]] + +# objective for optimization +initial_gate = qeye(2) +target_gate = gates.hadamard_transform() + +times = np.linspace(0, np.pi / 2, 250) +``` + +## Crab algorithm + + +```python +n_params = 6 # adjust in steps of 3 +control_params = { + "ctrl_x": {"guess": [1 for _ in range(n_params)], "bounds": [(-1, 1)] * n_params}, + "ctrl_y": {"guess": [1 for _ in range(n_params)], "bounds": [(-1, 1)] * n_params}, + "ctrl_z": {"guess": [1 for _ in range(n_params)], "bounds": [(-1, 1)] * n_params}, +} +alg_args = {"alg": "CRAB", "fid_err_targ": 0.01} + +res_crab = optimize_pulses( + objectives = Objective(initial_gate, H, target_gate), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "CRAB", + "fid_err_targ": 0.01 + }, +) + +print('Infidelity: ', res_crab.infidelity) + +plt.plot(times, res_crab.optimized_controls[0], 'b', label='optimized pulse sx') +plt.plot(times, res_crab.optimized_controls[1], 'g', label='optimized pulse sy') +plt.plot(times, res_crab.optimized_controls[2], 'r', label='optimized pulse sz') +plt.title('CRAB pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +initial_super = qt.to_super(initial_gate) +target_super = qt.to_super(target_gate) + +H_result = [Hd, + [Hc[0], np.array(res_crab.optimized_controls[0])], + [Hc[1], np.array(res_crab.optimized_controls[1])], + [Hc[2], np.array(res_crab.optimized_controls[2])]] +evolution = qt.mesolve(H_result, initial_super, times) + +plt.plot(times, [fidelity(gate, initial_super) for gate in evolution.states], label="Overlap with initial gate") +plt.plot(times, [fidelity(gate, target_super) for gate in evolution.states], label="Overlap with target gate") + +plt.title("CRAB performance") +plt.xlabel("Time") +plt.legend() +plt.show() +``` +## Validation + + +```python +assert res_crab.infidelity < 0.01 +assert np.allclose(fidelity(evolution.states[-1], target_super), 1 - res_crab.infidelity, atol=1e-3) +``` + + +```python +qt.about() +``` diff --git a/tests/interactive/CRAB_state_closed.md b/tests/interactive/CRAB_state_closed.md index 278ab21..5eff137 100644 --- a/tests/interactive/CRAB_state_closed.md +++ b/tests/interactive/CRAB_state_closed.md @@ -5,14 +5,14 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.17.1 + jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# CRAB algorithm for 2 level system +# CRAB algorithm for a closed system (state transfer) ```python import matplotlib.pyplot as plt @@ -87,7 +87,7 @@ plt.show() ```python assert res_crab.infidelity < 0.001 -assert np.abs(evolution.states[-1].overlap(target_state)) > 1-0.001 +assert np.allclose(np.abs(evolution.states[-1].overlap(target_state)), 1 - res_crab.infidelity, atol=1e-3) ``` ```python diff --git a/tests/interactive/CRAB_state_open.md b/tests/interactive/CRAB_state_open.md new file mode 100644 index 0000000..c26a9a9 --- /dev/null +++ b/tests/interactive/CRAB_state_open.md @@ -0,0 +1,110 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.17.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# CRAB algorithm for an open system (state transfer) + +```python +import matplotlib.pyplot as plt +import numpy as np +from qutip import basis, ket2dm, liouvillian, sigmam, Qobj +import qutip as qt +from qutip_qoc import Objective, optimize_pulses + +def fidelity(dm, target_dm): + """ + Fidelity used for density matrices in qutip-qtrl and qutip-qoc + """ + + diff = dm - target_dm + return 1 - np.real(diff.overlap(diff) / target_dm.norm()) / 2 +``` + +## Problem setup + + +```python +# Energy levels +E1, E2 = 1.0, 2.0 + +gamma = 0.1 # amplitude damping +c_ops = [np.sqrt(gamma) * sigmam()] + +Hd = Qobj(np.diag([E1, E2])) +Hd = liouvillian(H=Hd, c_ops=c_ops) +Hc = Qobj(np.array([ + [0, 1], + [1, 0] +])) +Hc = liouvillian(Hc) +H = [Hd, Hc] + +initial_state = ket2dm(basis(2, 0)) +target_state = ket2dm(basis(2, 1)) + +times = np.linspace(0, 2 * np.pi, 250) +``` + +## CRAB algorithm + + +```python +n_params = 6 # adjust in steps of 3 +control_params = { + "ctrl_x": {"guess": [1 for _ in range(n_params)], "bounds": [(-1, 1)] * n_params}, +} + +res_crab = optimize_pulses( + objectives = Objective(initial_state, H, target_state), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "CRAB", + "fid_err_targ": 0.02 + }, +) + +print('Infidelity: ', res_crab.infidelity) + +plt.plot(times, res_crab.optimized_controls[0], label='optimized pulse') +plt.title('CRAB pulse') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +H_result = [Hd, [Hc, res_crab.optimized_controls[0]]] +evolution = qt.mesolve(H_result, initial_state, times) + +plt.plot(times, [fidelity(dm, initial_state) for dm in evolution.states], label="Overlap with initial state") +plt.plot(times, [fidelity(dm, target_state) for dm in evolution.states], label="Overlap with target state") + +plt.title("CRAB performance") +plt.xlabel('Time') +plt.legend() +plt.show() +``` +## Validation + + +```python +assert res_crab.infidelity < 0.02 +assert np.allclose(fidelity(evolution.states[-1], target_state), 1 - res_crab.infidelity, atol=1e-3) +``` + + +```python +qt.about() +``` \ No newline at end of file diff --git a/tests/interactive/GOAT_gate_closed.md b/tests/interactive/GOAT_gate_closed.md index 4eafa97..9ffcc37 100644 --- a/tests/interactive/GOAT_gate_closed.md +++ b/tests/interactive/GOAT_gate_closed.md @@ -5,14 +5,14 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.17.1 + jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# GOAT algorithm for a closed system +# GOAT algorithm for a closed system (gate synthesis)| ```python import matplotlib.pyplot as plt @@ -187,7 +187,7 @@ plt.legend() plt.show() ``` -## Global optimization +### c) global optimization ```python res_goat_global = optimize_pulses( @@ -266,13 +266,13 @@ plt.show() ```python assert res_goat.infidelity < 0.001 -assert fidelity(evolution.states[-1], target_gate) > 1-0.001 +assert np.allclose(fidelity(evolution.states[-1], target_gate), 1 - res_goat.infidelity, atol=1e-3) assert res_goat_time.infidelity < 0.001 -assert fidelity(evolution_time.states[-1], target_gate) > 1-0.001 +assert np.allclose(fidelity(evolution_time.states[-1], target_gate), 1 - res_goat_time.infidelity, atol=1e-3) assert res_goat_global.infidelity < 0.001 -assert fidelity(evolution_global.states[-1], target_gate) > 1-0.001 +assert np.allclose(fidelity(evolution_global.states[-1], target_gate), 1 - res_goat_global.infidelity, atol=1e-3) ``` ```python diff --git a/tests/interactive/GOAT_gate_open.md b/tests/interactive/GOAT_gate_open.md new file mode 100644 index 0000000..1dca173 --- /dev/null +++ b/tests/interactive/GOAT_gate_open.md @@ -0,0 +1,289 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.17.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# GOAT algorithm for an open system (gate synthesis) + +```python +import matplotlib.pyplot as plt +import numpy as np +from qutip import gates, qeye, liouvillian, sigmam, sigmax, sigmay, sigmaz +import qutip as qt +from qutip_qoc import Objective, optimize_pulses + +def fidelity(gate_super, target_super): + gate_oper = qt.Qobj(gate_super.data) + target_oper = qt.Qobj(target_super.data) + + return np.abs(gate_oper.overlap(target_oper) / target_oper.norm()) +``` + +## Problem setup + + +```python +omega = 0.1 # energy splitting +gamma = 0.1 # amplitude damping +sx, sy, sz = sigmax(), sigmay(), sigmaz() +c_ops = [np.sqrt(gamma) * sigmam()] + +Hd = 1 / 2 * omega * sz +Hd = liouvillian(H=Hd, c_ops=c_ops) +Hc = [liouvillian(sx), liouvillian(sy), liouvillian(sz)] +H = [Hd, Hc[0], Hc[1], Hc[2]] + +# objective for optimization +initial_gate = qeye(2) +target_gate = gates.hadamard_transform() + +times = np.linspace(0, np.pi / 2, 250) +``` + +## Guess + + +```python +goat_guess = [1, 1] +guess_pulse = goat_guess[0] * np.sin(goat_guess[1] * times) + +initial_super = qt.to_super(initial_gate) +target_super = qt.to_super(target_gate) + +H_guess = [Hd] + [[hc, guess_pulse] for hc in Hc] +evolution_guess = qt.mesolve(H_guess, initial_super, times) + +plt.plot(times, [fidelity(gate, initial_super) for gate in evolution_guess.states], label="Overlap with initial gate") +plt.plot(times, [fidelity(gate, target_super) for gate in evolution_guess.states], label="Overlap with target gate") +plt.title("Guess performance") +plt.xlabel('Time') +plt.legend() +plt.show() +``` +## GOAT algorithm +### a) not optimized over time + +```python +# control function +def sin(t, c): + return c[0] * np.sin(c[1] * t) + +# derivatives +def grad_sin(t, c, idx): + if idx == 0: # w.r.t. c0 + return np.sin(c[1] * t) + if idx == 1: # w.r.t. c1 + return c[0] * np.cos(c[1] * t) * t + if idx == 2: # w.r.t. time + return c[0] * np.cos(c[1] * t) * c[1] + +H = [Hd] + [[hc, sin, {"grad": grad_sin}] for hc in Hc] +``` + + +```python +control_params = { + id: {"guess": goat_guess, "bounds": [(-1, 1), (0, 2 * np.pi)]} # c0 and c1 + for id in ['x', 'y', 'z'] +} + +# run the optimization +res_goat = optimize_pulses( + objectives = Objective(initial_gate, H, target_gate), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "GOAT", + "fid_err_targ": 0.001, + }, +) + +print('Infidelity: ', res_goat.infidelity) + +plt.plot(times, res_goat.optimized_controls[0], 'b', label='optimized pulse sx') +plt.plot(times, res_goat.optimized_controls[1], 'g', label='optimized pulse sy') +plt.plot(times, res_goat.optimized_controls[2], 'r', label='optimized pulse sz') +plt.title('GOAT pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +H_result = [Hd, + [Hc[0], np.array(res_goat.optimized_controls[0])], + [Hc[1], np.array(res_goat.optimized_controls[1])], + [Hc[2], np.array(res_goat.optimized_controls[2])]] +evolution = qt.mesolve(H_result, initial_super, times) + +plt.plot(times, [fidelity(gate, initial_super) for gate in evolution.states], label="Overlap with initial gate") +plt.plot(times, [fidelity(gate, target_super) for gate in evolution.states], label="Overlap with target gate") + +plt.title("GOAT performance") +plt.xlabel("Time") +plt.legend() +plt.show() +``` +### b) optimized over time + + +```python +# treats time as optimization variable +control_params["__time__"] = { + "guess": times[len(times) // 2], + "bounds": [times[0], times[-1]], +} + +# run the optimization +res_goat_time = optimize_pulses( + objectives = Objective(initial_gate, H, target_gate), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "GOAT", + "fid_err_targ": 0.001, + }, +) + +opt_time = res_goat_time.optimized_params[-1][0] +time_range = times < opt_time + +print('Infidelity: ', res_goat_time.infidelity) +print('Optimized time: ', opt_time) + +plt.plot(times, guess_pulse, 'k--', label='guess pulse sx, sy, sz') +plt.plot(times[time_range], np.array(res_goat_time.optimized_controls[0])[time_range], 'b', label='optimized pulse sx') +plt.plot(times[time_range], np.array(res_goat_time.optimized_controls[1])[time_range], 'g', label='optimized pulse sy') +plt.plot(times[time_range], np.array(res_goat_time.optimized_controls[2])[time_range], 'r', label='optimized pulse sz') +plt.title('GOAT pulses (time optimization)') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +times2 = times[time_range] +if opt_time not in times2: + times2 = np.append(times2, opt_time) + +H_result = qt.QobjEvo( + [Hd, [Hc[0], np.array(res_goat_time.optimized_controls[0])], + [Hc[1], np.array(res_goat_time.optimized_controls[1])], + [Hc[2], np.array(res_goat_time.optimized_controls[2])]], tlist=times) +evolution_time = qt.mesolve(H_result, initial_super, times2) + +plt.plot(times2, [fidelity(gate, initial_super) for gate in evolution_time.states], label="Overlap with initial gate") +plt.plot(times2, [fidelity(gate, target_super) for gate in evolution_time.states], label="Overlap with target gate") + +plt.title('GOAT (optimized over time) performance') +plt.xlabel("Time") +plt.legend() +plt.show() +``` +### c) global optimization + +```python +res_goat_global = optimize_pulses( + objectives = Objective(initial_gate, H, target_gate), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "GOAT", + "fid_err_targ": 0.001, + }, + optimizer_kwargs = { + "method": "basinhopping", + "max_iter": 100, + } +) + +global_time = res_goat_global.optimized_params[-1][0] +global_range = times < global_time + +print('Infidelity: ', res_goat_global.infidelity) +print('Optimized time: ', global_time) + +plt.plot(times, guess_pulse, 'k--', label='guess pulse sx, sy, sz') +plt.plot(times[global_range], np.array(res_goat_global.optimized_controls[0])[global_range], 'b', label='optimized pulse sx') +plt.plot(times[global_range], np.array(res_goat_global.optimized_controls[1])[global_range], 'g', label='optimized pulse sy') +plt.plot(times[global_range], np.array(res_goat_global.optimized_controls[2])[global_range], 'r', label='optimized pulse sz') +plt.title('GOAT pulses (global optimization)') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +times3 = times[global_range] +if global_time not in times3: + times3 = np.append(times3, global_time) + +H_result = qt.QobjEvo( + [Hd, [Hc[0], np.array(res_goat_global.optimized_controls[0])], + [Hc[1], np.array(res_goat_global.optimized_controls[1])], + [Hc[2], np.array(res_goat_global.optimized_controls[2])]], tlist=times) +evolution_global = qt.mesolve(H_result, initial_super, times3) + +plt.plot(times3, [fidelity(gate, initial_super) for gate in evolution_global.states], label="Overlap with initial gate") +plt.plot(times3, [fidelity(gate, target_super) for gate in evolution_global.states], label="Overlap with target gate") + +plt.title('GOAT (global optimization) performance') +plt.xlabel('Time') +plt.legend() +plt.show() +``` +## Comparison + + +```python +fig, axes = plt.subplots(1, 3, figsize=(18, 4)) # 1 row, 3 columns + +titles = ["GOAT sx pulses", "GOAT sy pulses", "GOAT sz pulses"] + +for i, ax in enumerate(axes): + ax.plot(times, guess_pulse, label='initial guess') + ax.plot(times, res_goat.optimized_controls[i], label='optimized pulse') + ax.plot(times[time_range], np.array(res_goat_time.optimized_controls[i])[time_range], label='optimized (over time) pulse') + ax.plot(times[global_range], np.array(res_goat_global.optimized_controls[i])[global_range], label='global optimized pulse') + ax.set_title(titles[i]) + ax.set_xlabel('Time') + ax.set_ylabel('Pulse amplitude') + ax.legend() + +plt.tight_layout() +plt.show() +``` +## Validation + + +```python +guess_fidelity = fidelity(evolution_guess.states[-1], target_super) + +# target fidelity not reached in part a), check that it is better than the guess +assert 1 - res_goat.infidelity >= guess_fidelity +assert np.allclose(fidelity(evolution.states[-1], target_super), 1 - res_goat.infidelity, atol=1e-3) + +# target fidelity not reached in part b), check that it is better than part a) +assert res_goat_time.infidelity <= res_goat.infidelity +assert np.allclose(fidelity(evolution_time.states[-1], target_super), 1 - res_goat_time.infidelity, atol=1e-3) + +assert res_goat_global.infidelity <= res_goat_time.infidelity +assert np.allclose(fidelity(evolution_global.states[-1], target_super), 1 - res_goat_global.infidelity, atol=1e-3) +``` + + +```python +qt.about() +``` \ No newline at end of file diff --git a/tests/interactive/GOAT_state_closed.md b/tests/interactive/GOAT_state_closed.md index bbce5f2..e0316b4 100644 --- a/tests/interactive/GOAT_state_closed.md +++ b/tests/interactive/GOAT_state_closed.md @@ -5,14 +5,14 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.17.1 + jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# GOAT algorithm for a 2 level system +# GOAT algorithm for a closed system (state transfer) ```python import matplotlib.pyplot as plt @@ -284,7 +284,7 @@ assert res_goat_time.infidelity < res_goat.infidelity assert np.allclose(np.abs(evolution_time.states[-1].overlap(target_state)), 1 - res_goat_time.infidelity, atol=1e-3) assert res_goat_global.infidelity < 0.001 -assert np.abs(evolution_global.states[-1].overlap(target_state)) > 1 - 0.001 +assert np.allclose(np.abs(evolution_global.states[-1].overlap(target_state)), 1 - res_goat_global.infidelity, atol=1e-3) ``` ```python diff --git a/tests/interactive/GOAT_state_open.md b/tests/interactive/GOAT_state_open.md new file mode 100644 index 0000000..e92e900 --- /dev/null +++ b/tests/interactive/GOAT_state_open.md @@ -0,0 +1,298 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.17.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# GOAT algorithm for an open system (state transfer) + +```python +import matplotlib.pyplot as plt +import numpy as np +from qutip import basis, ket2dm, liouvillian, sigmam, Qobj +import qutip as qt +from qutip_qoc import Objective, optimize_pulses + +def fidelity(dm, target_dm): + """ + Fidelity used for density matrices in qutip-qtrl and qutip-qoc + """ + + diff = dm - target_dm + return 1 - np.real(diff.overlap(diff) / target_dm.norm()) / 2 +``` + +## Problem setup + + +```python +# Energy levels +E1, E2 = 1.0, 2.0 + +gamma = 0.1 # amplitude damping +c_ops = [np.sqrt(gamma) * sigmam()] + +Hd = Qobj(np.diag([E1, E2])) +Hd = liouvillian(H=Hd, c_ops=c_ops) +Hc = Qobj(np.array([ + [0, 1], + [1, 0] +])) +Hc = liouvillian(Hc) +H = [Hd, Hc] + +initial_state = ket2dm(basis(2, 0)) +target_state = ket2dm(basis(2, 1)) + +times = np.linspace(0, 2 * np.pi, 250) +``` + +## Guess + + +```python +goat_guess = [1, 0.5] +guess_pulse = goat_guess[0] * np.sin(goat_guess[1] * times) + +H_guess = [Hd, [Hc, guess_pulse]] +evolution_guess = qt.mesolve(H_guess, initial_state, times) + +print('Fidelity: ', fidelity(evolution_guess.states[-1], target_state)) + +plt.plot(times, [fidelity(state, initial_state) for state in evolution_guess.states], label="Overlap with initial state") +plt.plot(times, [fidelity(state, target_state) for state in evolution_guess.states], label="Overlap with target state") +plt.title("Guess performance") +plt.xlabel("Time") +plt.legend() +plt.show() +``` + +## GOAT algorithm + + +```python +# control function +def sin(t, c): + return c[0] * np.sin(c[1] * t) + +# gradient +def grad_sin(t, c, idx): + if idx == 0: # w.r.t. c0 + return np.sin(c[1] * t) + if idx == 1: # w.r.t. c1 + return c[0] * np.cos(c[1] * t) * t + if idx == 2: # w.r.t. time + return c[0] * np.cos(c[1] * t) * c[1] + +H = [Hd] + [[Hc, sin, {"grad": grad_sin}]] +``` + +### a) not optimized over time + + +```python +control_params = { + "ctrl_x": {"guess": goat_guess, "bounds": [(-3, 3), (0, 2 * np.pi)]} # c0 and c1 +} + +# run the optimization +res_goat = optimize_pulses( + objectives = Objective(initial_state, H, target_state), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "GOAT", + "fid_err_targ": 0.01 + }, +) + +print('Infidelity: ', res_goat.infidelity) + +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_goat.optimized_controls[0], label='optimized pulse') +plt.title('GOAT pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +H_result = [Hd, [Hc, np.array(res_goat.optimized_controls[0])]] +evolution = qt.mesolve(H_result, initial_state, times) + +plt.plot(times, [fidelity(state, initial_state) for state in evolution.states], label="Overlap with initial state") +plt.plot(times, [fidelity(state, target_state) for state in evolution.states], label="Overlap with target state") + +plt.title("GOAT performance") +plt.xlabel("Time") +plt.legend() +plt.show() +``` +The desired fidelity is not reached. + + +### b) optimized over time + + +```python +# treats time as optimization variable +control_params["__time__"] = { + "guess": times[len(times) // 2], + "bounds": [times[0], times[-1]], +} + +# run the optimization +res_goat_time = optimize_pulses( + objectives = Objective(initial_state, H, target_state), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "GOAT", + "fid_err_targ": 0.01 + }, +) + +opt_time = res_goat_time.optimized_params[-1][0] +time_range = times < opt_time + +print('Infidelity: ', res_goat_time.infidelity) +print('Optimized time: ', opt_time) + +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_goat.optimized_controls[0], label='optimized pulse') +plt.plot(times[time_range], np.array(res_goat_time.optimized_controls[0])[time_range], label='optimized (over time) pulse') +plt.title('GOAT pulses (time optimization)') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +times2 = times[time_range] +if opt_time not in times2: + times2 = np.append(times2, opt_time) + +H_result = qt.QobjEvo([Hd, [Hc, np.array(res_goat_time.optimized_controls[0])]], tlist=times) +evolution_time = qt.mesolve(H_result, initial_state, times2) + +plt.plot(times2, [fidelity(state, initial_state) for state in evolution_time.states], label="Overlap with initial state") +plt.plot(times2, [fidelity(state, target_state) for state in evolution_time.states], label="Overlap with target state") + +plt.title('GOAT (optimized over time) performance') +plt.xlabel('Time') +plt.legend() +plt.show() +``` +The desired fidelity is still not reached. + + +### c) global optimization + +```python +res_goat_global = optimize_pulses( + objectives = Objective(initial_state, H, target_state), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "GOAT", + "fid_err_targ": 0.01 + }, + optimizer_kwargs={ + "method": "basinhopping", + "max_iter": 1000 + } +) + +global_time = res_goat_global.optimized_params[-1][0] +global_range = times < global_time + +print('Infidelity: ', res_goat_global.infidelity) +print('Optimized time: ', global_time) + +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_goat.optimized_controls[0], label='optimized pulse') +plt.plot(times[global_range], np.array(res_goat_global.optimized_controls[0])[global_range], label='global optimized pulse') +plt.title('GOAT pulses (global)') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +times3 = times[global_range] +if global_time not in times3: + times3 = np.append(times3, global_time) + +H_result = qt.QobjEvo([Hd, [Hc, np.array(res_goat_global.optimized_controls[0])]], tlist=times) +evolution_global = qt.mesolve(H_result, initial_state, times3) + +plt.plot(times3, [fidelity(state, initial_state) for state in evolution_global.states], label="Overlap with initial state") +plt.plot(times3, [fidelity(state, target_state) for state in evolution_global.states], label="Overlap with target state") + +plt.title('GOAT (global) performance') +plt.xlabel('Time') +plt.legend() +plt.show() +``` +## Comparison + +```python +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_goat.optimized_controls[0], label='optimized pulse') +plt.plot(times[time_range], np.array(res_goat_time.optimized_controls[0])[time_range], label='optimized (over time) pulse') +plt.plot(times[global_range], np.array(res_goat_global.optimized_controls[0])[global_range], label='global optimized pulse') +plt.title('GOAT pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` +```python +print('Guess Fidelity: ', fidelity(evolution_guess.states[-1], target_state)) +print('GOAT Fidelity: ', 1 - res_goat.infidelity) +print('Time Fidelity: ', 1 - res_goat_time.infidelity) +print('GLobal Fidelity: ', 1 - res_goat_global.infidelity) + +plt.plot(times, [qt.fidelity(state, target_state) for state in evolution_guess.states], 'k--', label="Guess") +plt.plot(times, [qt.fidelity(state, target_state) for state in evolution.states], label="GOAT") +plt.plot(times2, [qt.fidelity(state, target_state) for state in evolution_time.states], label="Time") +plt.plot(times3, [qt.fidelity(state, target_state) for state in evolution_global.states], label="Global") + +plt.title('GOAT Fidelities') +plt.xlabel('Time') +plt.legend() +plt.show() +``` + +## Validation + + +```python +guess_fidelity = fidelity(evolution_guess.states[-1], target_state) + +# target fidelity not reached in part a), check that it is better than the guess +assert 1 - res_goat.infidelity > guess_fidelity +assert np.allclose(fidelity(evolution.states[-1], target_state), 1 - res_goat.infidelity, atol=1e-3) + +# target fidelity not reached in part b), check that it is better than part a) +assert res_goat_time.infidelity < res_goat.infidelity +assert np.allclose(fidelity(evolution_time.states[-1], target_state), 1 - res_goat_time.infidelity, atol=1e-3) + +assert res_goat_global.infidelity < 0.01 +assert np.allclose(fidelity(evolution_global.states[-1], target_state), 1 - res_goat_global.infidelity, atol=1e-3) +``` + +```python +qt.about() +``` \ No newline at end of file diff --git a/tests/interactive/GRAPE_gate_closed.md b/tests/interactive/GRAPE_gate_closed.md index bb9b145..43d5d5c 100644 --- a/tests/interactive/GRAPE_gate_closed.md +++ b/tests/interactive/GRAPE_gate_closed.md @@ -5,14 +5,14 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.17.1 + jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# GRAPE algorithm for a closed system +# GRAPE algorithm for a closed system (gate synthesis) ```python import matplotlib.pyplot as plt @@ -69,9 +69,9 @@ plt.show() ```python control_params = { - "ctrl_x": {"guess": np.sin(times), "bounds": [-1, 1]}, - "ctrl_y": {"guess": np.cos(times), "bounds": [-1, 1]}, - "ctrl_z": {"guess": np.tanh(times), "bounds": [-1, 1]}, + "ctrl_x": {"guess": guess_pulse_x, "bounds": [-1, 1]}, + "ctrl_y": {"guess": guess_pulse_y, "bounds": [-1, 1]}, + "ctrl_z": {"guess": guess_pulse_z, "bounds": [-1, 1]}, } res_grape = optimize_pulses( @@ -100,7 +100,10 @@ plt.show() ``` ```python -H_result = [Hd, [Hc[0], res_grape.optimized_controls[0]], [Hc[1], res_grape.optimized_controls[1]], [Hc[2], res_grape.optimized_controls[2]]] +H_result = [Hd, + [Hc[0], np.array(res_grape.optimized_controls[0])], + [Hc[1], np.array(res_grape.optimized_controls[1])], + [Hc[2], np.array(res_grape.optimized_controls[2])]] evolution = qt.sesolve(H_result, initial_gate, times) plt.plot(times, [fidelity(gate, initial_gate) for gate in evolution.states], label="Overlap with initial gate") @@ -116,7 +119,7 @@ plt.show() ```python assert res_grape.infidelity < 0.001 -assert fidelity(evolution.states[-1], target_gate) > 1-0.001 +assert np.allclose(fidelity(evolution.states[-1], target_gate), 1 - res_grape.infidelity, atol=1e-3) ``` ```python diff --git a/tests/interactive/GRAPE_gate_open.md b/tests/interactive/GRAPE_gate_open.md new file mode 100644 index 0000000..50e3807 --- /dev/null +++ b/tests/interactive/GRAPE_gate_open.md @@ -0,0 +1,136 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.17.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# GRAPE algorithm for an open system (gate synthesis) + +```python +import matplotlib.pyplot as plt +import numpy as np +from qutip import gates, qeye, liouvillian, sigmam, sigmax, sigmay, sigmaz +import qutip as qt +from qutip_qoc import Objective, optimize_pulses + +def fidelity(gate_super, target_super): + gate_oper = qt.Qobj(gate_super.data) + target_oper = qt.Qobj(target_super.data) + + return np.abs(gate_oper.overlap(target_oper) / target_oper.norm()) +``` + +## Problem setup + + +```python +omega = 0.1 # energy splitting +gamma = 0.1 # amplitude damping +sx, sy, sz = sigmax(), sigmay(), sigmaz() +c_ops = [np.sqrt(gamma) * sigmam()] + +Hd = 1 / 2 * omega * sz +Hd = liouvillian(H=Hd, c_ops=c_ops) +Hc = [liouvillian(sx), liouvillian(sy), liouvillian(sz)] +H = [Hd, Hc[0], Hc[1], Hc[2]] + +# objective for optimization +initial_gate = qeye(2) +target_gate = gates.hadamard_transform() + +times = np.linspace(0, np.pi / 2, 250) +``` + +## Guess + + +```python +guess_pulse_x = np.sin(times) +guess_pulse_y = np.cos(times) +guess_pulse_z = np.tanh(times) + +initial_super = qt.to_super(initial_gate) +target_super = qt.to_super(target_gate) + +H_guess = [Hd, [Hc[0], guess_pulse_x], [Hc[1], guess_pulse_y], [Hc[2], guess_pulse_z]] +evolution_guess = qt.mesolve(H_guess, initial_super, times) + +print('Fidelity: ', fidelity(evolution_guess.states[-1], target_super)) + +plt.plot(times, [fidelity(gate, initial_super) for gate in evolution_guess.states], label="Overlap with initial gate") +plt.plot(times, [fidelity(gate, target_super) for gate in evolution_guess.states], label="Overlap with target gate") +plt.title("Guess performance") +plt.xlabel('Time') +plt.legend() +plt.show() +``` +## Grape algorithm + + +```python +control_params = { + "ctrl_x": {"guess": guess_pulse_x, "bounds": [-1, 1]}, + "ctrl_y": {"guess": guess_pulse_y, "bounds": [-1, 1]}, + "ctrl_z": {"guess": guess_pulse_z, "bounds": [-1, 1]}, +} + +res_grape = optimize_pulses( + objectives = Objective(initial_gate, H, target_gate), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "GRAPE", + "fid_err_targ": 0.01 + }, +) + +print('Infidelity: ', res_grape.infidelity) + +plt.plot(times, guess_pulse_x, 'b--', label='guess pulse sx') +plt.plot(times, res_grape.optimized_controls[0], 'b', label='optimized pulse sx') +plt.plot(times, guess_pulse_y, 'g--', label='guess pulse sy') +plt.plot(times, res_grape.optimized_controls[1], 'g', label='optimized pulse sy') +plt.plot(times, guess_pulse_z, 'r--', label='guess pulse sz') +plt.plot(times, res_grape.optimized_controls[2], 'r', label='optimized pulse sz') +plt.title('GRAPE pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +H_result = [Hd, + [Hc[0], np.array(res_grape.optimized_controls[0])], + [Hc[1], np.array(res_grape.optimized_controls[1])], + [Hc[2], np.array(res_grape.optimized_controls[2])]] +evolution = qt.mesolve(H_result, initial_super, times) + +plt.plot(times, [fidelity(gate, initial_super) for gate in evolution.states], label="Overlap with initial gate") +plt.plot(times, [fidelity(gate, target_super) for gate in evolution.states], label="Overlap with target gate") + +plt.title('GRAPE performance') +plt.xlabel('Time') +plt.legend() +plt.show() +``` +## Validation + + +```python +assert res_grape.infidelity < 0.01 +assert np.allclose(fidelity(evolution.states[-1], target_super), 1 - res_grape.infidelity, atol=1e-3) +``` + + +```python +qt.about() +``` \ No newline at end of file diff --git a/tests/interactive/GRAPE_state_closed.md b/tests/interactive/GRAPE_state_closed.md index 2681df7..7893a17 100644 --- a/tests/interactive/GRAPE_state_closed.md +++ b/tests/interactive/GRAPE_state_closed.md @@ -5,14 +5,14 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.17.1 + jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# GRAPE algorithm for 2 level system +# GRAPE algorithm for a closed system (state transfer) ```python import matplotlib.pyplot as plt @@ -106,7 +106,7 @@ plt.show() ```python assert res_grape.infidelity < 0.01 -assert np.abs(evolution.states[-1].overlap(target_state)) > 1-0.01 +assert np.allclose(np.abs(evolution.states[-1].overlap(target_state)), 1 - res_grape.infidelity, atol=1e-3) ``` ```python diff --git a/tests/interactive/GRAPE_state_open.md b/tests/interactive/GRAPE_state_open.md new file mode 100644 index 0000000..c9b05dd --- /dev/null +++ b/tests/interactive/GRAPE_state_open.md @@ -0,0 +1,131 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.17.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# GRAPE algorithm for an open system (state transfer) + +```python +import matplotlib.pyplot as plt +import numpy as np +from qutip import basis, ket2dm, liouvillian, sigmam, Qobj +import qutip as qt +from qutip_qoc import Objective, optimize_pulses + +def fidelity(dm, target_dm): + """ + Fidelity used for density matrices in qutip-qtrl and qutip-qoc + """ + + diff = dm - target_dm + return 1 - np.real(diff.overlap(diff) / target_dm.norm()) / 2 +``` + +## Problem setup + + +```python +# Energy levels +E1, E2 = 1.0, 2.0 + +gamma = 0.1 # amplitude damping +c_ops = [np.sqrt(gamma) * sigmam()] + +Hd = Qobj(np.diag([E1, E2])) +Hd = liouvillian(H=Hd, c_ops=c_ops) +Hc = Qobj(np.array([ + [0, 1], + [1, 0] +])) +Hc = liouvillian(Hc) +H = [Hd, Hc] + +initial_state = ket2dm(basis(2, 0)) +target_state = ket2dm(basis(2, 1)) + +times = np.linspace(0, 2 * np.pi, 250) +``` + +## Guess + + +```python +guess_pulse = np.sin(times) + +H_guess = [Hd, [Hc, guess_pulse]] +evolution_guess = qt.mesolve(H_guess, initial_state, times) + +print('Fidelity: ', fidelity(evolution_guess.states[-1], target_state)) + +plt.plot(times, [fidelity(state, initial_state) for state in evolution_guess.states], label="Overlap with initial state") +plt.plot(times, [fidelity(state, target_state) for state in evolution_guess.states], label="Overlap with target state") +plt.title("Guess performance") +plt.xlabel("Time") +plt.legend() +plt.show() +``` + +## GRAPE algorithm + + +```python +control_params = { + "ctrl_x": {"guess": guess_pulse, "bounds": [-1, 1]}, # Control pulse for Hc1 +} + +res_grape = optimize_pulses( + objectives = Objective(initial_state, H, target_state), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "GRAPE", + "fid_err_targ": 0.01, + }, +) + +print('Infidelity: ', res_grape.infidelity) + +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_grape.optimized_controls[0], label='optimized pulse') +plt.title('GRAPE pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +control_step_function = qt.coefficient(res_grape.optimized_controls[0], tlist=times, order=0) +H_result = [Hd, [Hc, control_step_function]] + +fine_times = np.linspace(0, times[-1], len(times) * 25) +evolution = qt.mesolve(H_result, initial_state, fine_times) + +plt.plot(fine_times, [fidelity(state, initial_state) for state in evolution.states], label="Overlap with initial state") +plt.plot(fine_times, [fidelity(state, target_state) for state in evolution.states], label="Overlap with target state") +plt.title("GRAPE performance") +plt.xlabel('Time') +plt.legend() +plt.show() +``` +## Validation + + +```python +assert res_grape.infidelity < 0.01 +assert np.allclose(fidelity(evolution.states[-1], target_state), 1 - res_grape.infidelity, atol=1e-3) +``` + + +```python +qt.about() +``` diff --git a/tests/interactive/JOPT_gate_closed.md b/tests/interactive/JOPT_gate_closed.md index 1f3370c..dd3b335 100644 --- a/tests/interactive/JOPT_gate_closed.md +++ b/tests/interactive/JOPT_gate_closed.md @@ -5,14 +5,14 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.17.1 + jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# JOPT algorithm for a closed system +# JOPT algorithm for a closed system (gate synthesis) ```python import matplotlib.pyplot as plt @@ -190,7 +190,7 @@ plt.legend() plt.show() ``` -## Global optimization +### c) global optimization ```python res_jopt_global = optimize_pulses( @@ -269,13 +269,13 @@ plt.show() ```python assert res_jopt.infidelity < 0.001 -assert fidelity(evolution.states[-1], target_gate) > 1-0.001 +assert np.allclose(fidelity(evolution.states[-1], target_gate), 1 - res_jopt.infidelity, atol=1e-3) assert res_jopt_time.infidelity < 0.001 -assert fidelity(evolution_time.states[-1], target_gate) > 1-0.001 +assert np.allclose(fidelity(evolution_time.states[-1], target_gate), 1 - res_jopt_time.infidelity, atol=1e-3) assert res_jopt_global.infidelity < 0.001 -assert fidelity(evolution_global.states[-1], target_gate) > 1-0.001 +assert np.allclose(fidelity(evolution_global.states[-1], target_gate), 1 - res_jopt_global.infidelity, atol=1e-3) ``` ```python diff --git a/tests/interactive/JOPT_gate_open.md b/tests/interactive/JOPT_gate_open.md new file mode 100644 index 0000000..4a886e1 --- /dev/null +++ b/tests/interactive/JOPT_gate_open.md @@ -0,0 +1,298 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.17.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# JOPT algorithm for an open system (gate synthesis) + +```python +import matplotlib.pyplot as plt +import numpy as np +from qutip import gates, qeye, liouvillian, sigmam, sigmax, sigmay, sigmaz +import qutip as qt +from qutip_qoc import Objective, optimize_pulses + +try: + from jax import jit, numpy +except ImportError: # JAX not available, skip test + import pytest + pytest.skip("JAX not available") + +def fidelity(gate_super, target_super): + gate_oper = qt.Qobj(gate_super.data) + target_oper = qt.Qobj(target_super.data) + + return np.abs(gate_oper.overlap(target_oper) / target_oper.norm()) +``` + +## Problem setup + + +```python +omega = 0.1 # energy splitting +gamma = 0.1 # amplitude damping +sx, sy, sz = sigmax(), sigmay(), sigmaz() +c_ops = [np.sqrt(gamma) * sigmam()] + +Hd = 1 / 2 * omega * sz +Hd = liouvillian(H=Hd, c_ops=c_ops) +Hc = [liouvillian(sx), liouvillian(sy), liouvillian(sz)] +H = [Hd, Hc[0], Hc[1], Hc[2]] + +# objective for optimization +initial_gate = qeye(2) +target_gate = gates.hadamard_transform() + +times = np.linspace(0, np.pi / 2, 250) +``` + +## Guess + + +```python +jopt_guess = [1, 1] +guess_pulse = jopt_guess[0] * np.sin(jopt_guess[1] * times) + +initial_super = qt.to_super(initial_gate) +target_super = qt.to_super(target_gate) + +H_guess = [Hd] + [[hc, guess_pulse] for hc in Hc] +evolution_guess = qt.mesolve(H_guess, initial_super, times) + +plt.plot(times, [fidelity(gate, initial_super) for gate in evolution_guess.states], label="Overlap with initial gate") +plt.plot(times, [fidelity(gate, target_super) for gate in evolution_guess.states], label="Overlap with target gate") +plt.title("Guess performance") +plt.xlabel('Time') +plt.legend() +plt.show() +``` + + +## JOPT algorithm + + +```python +@jit +def sin_x(t, c, **kwargs): + return c[0] * numpy.sin(c[1] * t) + +H = [Hd] + [[hc, sin_x, {"grad": sin_x}] for hc in Hc] +``` + + +### a) not optimized over time + + +```python +control_params = { + id: {"guess": jopt_guess, "bounds": [(-1, 1), (0, 2 * np.pi)]} # c0 and c1 + for id in ['x', 'y', 'z'] +} + +# run the optimization +res_jopt = optimize_pulses( + objectives = Objective(initial_gate, H, target_gate), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "JOPT", + "fid_err_targ": 0.001, + }, +) + +print('Infidelity: ', res_jopt.infidelity) + +plt.plot(times, res_jopt.optimized_controls[0], 'b', label='optimized pulse sx') +plt.plot(times, res_jopt.optimized_controls[1], 'g', label='optimized pulse sy') +plt.plot(times, res_jopt.optimized_controls[2], 'r', label='optimized pulse sz') +plt.title('JOPT pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + + +```python +H_result = [Hd, + [Hc[0], np.array(res_jopt.optimized_controls[0])], + [Hc[1], np.array(res_jopt.optimized_controls[1])], + [Hc[2], np.array(res_jopt.optimized_controls[2])]] +evolution = qt.mesolve(H_result, initial_super, times) + +plt.plot(times, [fidelity(gate, initial_super) for gate in evolution.states], label="Overlap with initial gate") +plt.plot(times, [fidelity(gate, target_super) for gate in evolution.states], label="Overlap with target gate") + +plt.title("JOPT performance") +plt.xlabel("Time") +plt.legend() +plt.show() +``` + +### b) optimized over time + + +```python +# treats time as optimization variable +control_params["__time__"] = { + "guess": times[len(times) // 2], + "bounds": [times[0], times[-1]], +} + +# run the optimization +res_jopt_time = optimize_pulses( + objectives = Objective(initial_gate, H, target_gate), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "JOPT", + "fid_err_targ": 0.001, + }, +) + +opt_time = res_jopt_time.optimized_params[-1][0] +time_range = times < opt_time + +print('Infidelity: ', res_jopt_time.infidelity) +print('Optimized time: ', opt_time) + +plt.plot(times, guess_pulse, 'k--', label='guess pulse sx, sy, sz') +plt.plot(times[time_range], np.array(res_jopt_time.optimized_controls[0])[time_range], 'b', label='optimized pulse sx') +plt.plot(times[time_range], np.array(res_jopt_time.optimized_controls[1])[time_range], 'g', label='optimized pulse sy') +plt.plot(times[time_range], np.array(res_jopt_time.optimized_controls[2])[time_range], 'r', label='optimized pulse sz') +plt.title('JOPT pulses (time optimization)') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + + +```python +times2 = times[time_range] +if opt_time not in times2: + times2 = np.append(times2, opt_time) + +H_result = qt.QobjEvo( + [Hd, [Hc[0], np.array(res_jopt_time.optimized_controls[0])], + [Hc[1], np.array(res_jopt_time.optimized_controls[1])], + [Hc[2], np.array(res_jopt_time.optimized_controls[2])]], tlist=times) +evolution_time = qt.mesolve(H_result, initial_super, times2) + +plt.plot(times2, [fidelity(gate, initial_super) for gate in evolution_time.states], label="Overlap with initial gate") +plt.plot(times2, [fidelity(gate, target_super) for gate in evolution_time.states], label="Overlap with target gate") + +plt.title('JOPT (optimized over time) performance') +plt.xlabel("Time") +plt.legend() +plt.show() +``` + +### c) global optimization + +```python +res_jopt_global = optimize_pulses( + objectives = Objective(initial_gate, H, target_gate), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "JOPT", + "fid_err_targ": 0.001, + }, + optimizer_kwargs = { + "method": "basinhopping", + "max_iter": 100, + } +) + +global_time = res_jopt_global.optimized_params[-1][0] +global_range = times < global_time + +print('Infidelity: ', res_jopt_global.infidelity) +print('Optimized time: ', global_time) + +plt.plot(times, guess_pulse, 'k--', label='guess pulse sx, sy, sz') +plt.plot(times[global_range], np.array(res_jopt_global.optimized_controls[0])[global_range], 'b', label='optimized pulse sx') +plt.plot(times[global_range], np.array(res_jopt_global.optimized_controls[1])[global_range], 'g', label='optimized pulse sy') +plt.plot(times[global_range], np.array(res_jopt_global.optimized_controls[2])[global_range], 'r', label='optimized pulse sz') +plt.title('JOPT pulses (global optimization)') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + + +```python +times3 = times[global_range] +if global_time not in times3: + times3 = np.append(times3, global_time) + +H_result = qt.QobjEvo( + [Hd, [Hc[0], np.array(res_jopt_global.optimized_controls[0])], + [Hc[1], np.array(res_jopt_global.optimized_controls[1])], + [Hc[2], np.array(res_jopt_global.optimized_controls[2])]], tlist=times) +evolution_global = qt.mesolve(H_result, initial_super, times3) + +plt.plot(times3, [fidelity(gate, initial_super) for gate in evolution_global.states], label="Overlap with initial gate") +plt.plot(times3, [fidelity(gate, target_super) for gate in evolution_global.states], label="Overlap with target gate") + +plt.title('JOPT (global optimization) performance') +plt.xlabel('Time') +plt.legend() +plt.show() +``` + +## Comparison + + +```python +fig, axes = plt.subplots(1, 3, figsize=(18, 4)) # 1 row, 3 columns + +titles = ["JOPT sx pulses", "JOPT sy pulses", "JOPT sz pulses"] + +for i, ax in enumerate(axes): + ax.plot(times, guess_pulse, label='initial guess') + ax.plot(times, res_jopt.optimized_controls[i], label='optimized pulse') + ax.plot(times[time_range], np.array(res_jopt_time.optimized_controls[i])[time_range], label='optimized (over time) pulse') + ax.plot(times[global_range], np.array(res_jopt_global.optimized_controls[i])[global_range], label='global optimized pulse') + ax.set_title(titles[i]) + ax.set_xlabel('Time') + ax.set_ylabel('Pulse amplitude') + ax.legend() + +plt.tight_layout() +plt.show() +``` + +## Validation + + +```python +guess_fidelity = fidelity(evolution_guess.states[-1], target_super) + +# target fidelity not reached in part a), check that it is better than the guess +assert 1 - res_jopt.infidelity >= guess_fidelity +assert np.allclose(fidelity(evolution.states[-1], target_super), 1 - res_jopt.infidelity, atol=1e-3) + +# target fidelity not reached in part b), check that it is better than part a) +assert res_jopt_time.infidelity <= res_jopt.infidelity +assert np.allclose(fidelity(evolution_time.states[-1], target_super), 1 - res_jopt_time.infidelity, atol=1e-3) + +assert res_jopt_global.infidelity <= res_jopt_time.infidelity +assert np.allclose(fidelity(evolution_global.states[-1], target_super), 1 - res_jopt_global.infidelity, atol=1e-3) +``` + + +```python +qt.about() +``` diff --git a/tests/interactive/JOPT_state_closed.md b/tests/interactive/JOPT_state_closed.md index 11547cb..b440642 100644 --- a/tests/interactive/JOPT_state_closed.md +++ b/tests/interactive/JOPT_state_closed.md @@ -5,15 +5,14 @@ jupyter: extension: .md format_name: markdown format_version: '1.3' - jupytext_version: 1.17.1 + jupytext_version: 1.17.2 kernelspec: display_name: Python 3 (ipykernel) language: python name: python3 --- -# JOPT algorithm for a 2 level system - +# JOPT algorithm for a closed system (state transfer) ```python import matplotlib.pyplot as plt @@ -288,7 +287,7 @@ assert res_jopt_time.infidelity < res_jopt.infidelity assert np.allclose(np.abs(evolution_time.states[-1].overlap(target_state)), 1 - res_jopt_time.infidelity, atol=1e-3) assert res_jopt_global.infidelity < 0.001 -assert np.abs(evolution_global.states[-1].overlap(target_state)) > 1 - 0.001 +assert np.allclose(np.abs(evolution_global.states[-1].overlap(target_state)), 1 - res_jopt_global.infidelity, atol=1e-3) ``` ```python diff --git a/tests/interactive/JOPT_state_open.md b/tests/interactive/JOPT_state_open.md new file mode 100644 index 0000000..9392cbe --- /dev/null +++ b/tests/interactive/JOPT_state_open.md @@ -0,0 +1,291 @@ +--- +jupyter: + jupytext: + text_representation: + extension: .md + format_name: markdown + format_version: '1.3' + jupytext_version: 1.17.2 + kernelspec: + display_name: Python 3 (ipykernel) + language: python + name: python3 +--- + +# JOPT algorithm for an open system (state transfer) + +```python +import matplotlib.pyplot as plt +import numpy as np +from qutip import basis, ket2dm, liouvillian, sigmam, Qobj +import qutip as qt +from qutip_qoc import Objective, optimize_pulses + +try: + from jax import jit, numpy +except ImportError: # JAX not available, skip test + import pytest + pytest.skip("JAX not available") + +def fidelity(dm, target_dm): + """ + Fidelity used for density matrices in qutip-qtrl and qutip-qoc + """ + + diff = dm - target_dm + return 1 - np.real(diff.overlap(diff) / target_dm.norm()) / 2 +``` + +## Problem setup + + +```python +# Energy levels +E1, E2 = 1.0, 2.0 + +gamma = 0.1 # amplitude damping +c_ops = [np.sqrt(gamma) * sigmam()] + +Hd = Qobj(np.diag([E1, E2])) +Hd = liouvillian(H=Hd, c_ops=c_ops) +Hc = Qobj(np.array([ + [0, 1], + [1, 0] +])) +Hc = liouvillian(Hc) +H = [Hd, Hc] + +initial_state = ket2dm(basis(2, 0)) +target_state = ket2dm(basis(2, 1)) + +times = np.linspace(0, 2 * np.pi, 250) +``` + +## Guess + + +```python +jopt_guess = [1, 0.5] +guess_pulse = jopt_guess[0] * np.sin(jopt_guess[1] * times) + +H_guess = [Hd, [Hc, guess_pulse]] +evolution_guess = qt.mesolve(H_guess, initial_state, times) + +print('Fidelity: ', fidelity(evolution_guess.states[-1], target_state)) + +plt.plot(times, [fidelity(state, initial_state) for state in evolution_guess.states], label="Overlap with initial state") +plt.plot(times, [fidelity(state, target_state) for state in evolution_guess.states], label="Overlap with target state") +plt.title("Guess performance") +plt.xlabel("Time") +plt.legend() +plt.show() +``` + +## JOPT algorithm + + +```python +@jit +def sin_x(t, c, **kwargs): + return c[0] * numpy.sin(c[1] * t) + +H = [Hd, [Hc, sin_x]] +``` + +### a) not optimized over time + + +```python +control_params = { + "ctrl_x": {"guess": jopt_guess, "bounds": [(-3, 3), (0, 2 * np.pi)]} # c0 and c1 +} + +# run the optimization +res_jopt = optimize_pulses( + objectives = Objective(initial_state, H, target_state), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "JOPT", + "fid_err_targ": 0.01 + }, +) + +print('Infidelity: ', res_jopt.infidelity) + +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_jopt.optimized_controls[0], label='optimized pulse') +plt.title('JOPT pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +H_result = [Hd, [Hc, np.array(res_jopt.optimized_controls[0])]] +evolution = qt.mesolve(H_result, initial_state, times) + +plt.plot(times, [fidelity(state, initial_state) for state in evolution.states], label="Overlap with initial state") +plt.plot(times, [fidelity(state, target_state) for state in evolution.states], label="Overlap with target state") + +plt.title("JOPT performance") +plt.xlabel("Time") +plt.legend() +plt.show() +``` +### b) optimized over time + + +```python +# treats time as optimization variable +control_params["__time__"] = { + "guess": times[len(times) // 2], + "bounds": [times[0], times[-1]], +} + +# run the optimization +res_jopt_time = optimize_pulses( + objectives = Objective(initial_state, H, target_state), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "JOPT", + "fid_err_targ": 0.01 + }, +) + +opt_time = res_jopt_time.optimized_params[-1][0] +time_range = times < opt_time + +print('Infidelity: ', res_jopt_time.infidelity) +print('Optimized time: ', opt_time) + +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_jopt.optimized_controls[0], label='optimized pulse') +plt.plot(times[time_range], np.array(res_jopt_time.optimized_controls[0])[time_range], label='optimized (over time) pulse') +plt.title('JOPT pulses (time optimization)') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +times2 = times[time_range] +if opt_time not in times2: + times2 = np.append(times2, opt_time) + +H_result = qt.QobjEvo([Hd, [Hc, np.array(res_jopt_time.optimized_controls[0])]], tlist=times) +evolution_time = qt.mesolve(H_result, initial_state, times2) + +plt.plot(times2, [fidelity(state, initial_state) for state in evolution_time.states], label="Overlap with initial state") +plt.plot(times2, [fidelity(state, target_state) for state in evolution_time.states], label="Overlap with target state") + +plt.title('JOPT (optimized over time) performance') +plt.xlabel('Time') +plt.legend() +plt.show() +``` +### c) global optimization + +```python +res_jopt_global = optimize_pulses( + objectives = Objective(initial_state, H, target_state), + control_parameters = control_params, + tlist = times, + algorithm_kwargs = { + "alg": "JOPT", + "fid_err_targ": 0.01 + }, + optimizer_kwargs={ + "method": "basinhopping", + "max_iter": 1000 + } +) + +global_time = res_jopt_global.optimized_params[-1][0] +global_range = times < global_time + +print('Infidelity: ', res_jopt_global.infidelity) +print('Optimized time: ', global_time) + +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_jopt.optimized_controls[0], label='optimized pulse') +plt.plot(times[global_range], np.array(res_jopt_global.optimized_controls[0])[global_range], label='global optimized pulse') +plt.title('JOPT pulses (global)') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` + +```python +times3 = times[global_range] +if global_time not in times3: + times3 = np.append(times3, global_time) + +H_result = qt.QobjEvo([Hd, [Hc, np.array(res_jopt_global.optimized_controls[0])]], tlist=times) +evolution_global = qt.mesolve(H_result, initial_state, times3) + +plt.plot(times3, [fidelity(state, initial_state) for state in evolution_global.states], label="Overlap with initial state") +plt.plot(times3, [fidelity(state, target_state) for state in evolution_global.states], label="Overlap with target state") + +plt.title('JOPT (global) performance') +plt.xlabel('Time') +plt.legend() +plt.show() +``` +## Comparison + + +```python +plt.plot(times, guess_pulse, label='initial guess') +plt.plot(times, res_jopt.optimized_controls[0], label='optimized pulse') +plt.plot(times[time_range], np.array(res_jopt_time.optimized_controls[0])[time_range], label='optimized (over time) pulse') +plt.plot(times[global_range], np.array(res_jopt_global.optimized_controls[0])[global_range], label='global optimized pulse') +plt.title('JOPT pulses') +plt.xlabel('Time') +plt.ylabel('Pulse amplitude') +plt.legend() +plt.show() +``` +```python +print('Guess Fidelity: ', fidelity(evolution_guess.states[-1], target_state)) +print('JOPT Fidelity: ', 1 - res_jopt.infidelity) +print('Time Fidelity: ', 1 - res_jopt_time.infidelity) +print('GLobal Fidelity: ', 1 - res_jopt_global.infidelity) + +plt.plot(times, [qt.fidelity(state, target_state) for state in evolution_guess.states], 'k--', label="Guess") +plt.plot(times, [qt.fidelity(state, target_state) for state in evolution.states], label="GOAT") +plt.plot(times2, [qt.fidelity(state, target_state) for state in evolution_time.states], label="Time") +plt.plot(times3, [qt.fidelity(state, target_state) for state in evolution_global.states], label="Global") + +plt.title('Fidelities') +plt.xlabel('Time') +plt.legend() +plt.show() +``` + +## Validation + + +```python +guess_fidelity = fidelity(evolution_guess.states[-1], target_state) + +# target fidelity not reached in part a), check that it is better than the guess +assert 1 - res_jopt.infidelity > guess_fidelity +assert np.allclose(fidelity(evolution.states[-1], target_state), 1 - res_jopt.infidelity, atol=1e-3) + +# target fidelity not reached in part b), check that it is better than part a) +assert res_jopt_time.infidelity < res_jopt.infidelity +assert np.allclose(fidelity(evolution_time.states[-1], target_state), 1 - res_jopt_time.infidelity, atol=1e-3) + +assert res_jopt_global.infidelity < 0.01 +assert np.allclose(fidelity(evolution_global.states[-1], target_state), 1 - res_jopt_global.infidelity, atol=1e-3) +``` + + +```python +qt.about() +``` \ No newline at end of file