{
"cells": [
{
"cell_type": "markdown",
"metadata": {},
"source": [
"# How to optimize controls using user-defined basis functions\n",
"**Create optimized controls using arbitrary basis functions**\n",
"\n",
"Boulder Opal exposes a highly-flexible [optimization engine](https://docs.q-ctrl.com/references/qctrl/Functions/calculate_optimization.html) for general-purpose gradient-based optimization.\n",
"The pulses can be described in terms of optimizable linear combinations from a set of user-defined basis functions, which can greatly reduce the dimensionality of the optimization search space.\n",
"In this notebook we will use [Lorentz functions](https://en.wikipedia.org/wiki/Cauchy_distribution), although the same technique has also seen success with other bases, for example [Slepian functions](https://journals.aps.org/pra/abstract/10.1103/PhysRevA.97.062346).\n",
"You can also read the related user guides showing how to find optimal pulses using a [Fourier basis](https://docs.q-ctrl.com/boulder-opal/user-guides/how-to-optimize-controls-using-a-fourier-basis) or a\n",
"[Hann series basis](https://docs.q-ctrl.com/boulder-opal/user-guides/how-to-optimize-controls-using-a-hann-series-basis)."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Summary workflow\n",
"\n",
"### 1. Define basis function for signal composition in the graph\n",
"\n",
"The Boulder Opal optimization engine provides a [library of operations](https://docs.q-ctrl.com/references/qctrl/Graphs.html#graph-operations) for creating optimizable signals in a custom basis.\n",
"In particular, you can create an analytical pulse by creating a node with the `graph.identity_stf` operation, which returns an STF representing the function $f(t) = t$, and then apply [mathematical and arithmetic operations](https://docs.q-ctrl.com/references/qctrl/Graphs.html#mathematical-functions) to it in order to generate more complex functions.\n",
"\n",
"In the example below, the function `custom_optimizable_superposition` defines the composition of the pulse in terms of Lorentz functions.\n",
"\n",
"### 2. Execute graph-based optimization\n",
"\n",
"With the graph object created, you can run an optimization using the `qctrl.functions.calculate_optimization` function, providing the graph, the cost node, and the desired output nodes.\n",
"The function returns the results of the optimization."
]
},
{
"cell_type": "markdown",
"metadata": {},
"source": [
"## Example: Robust optimization on a qutrit using a function superposition\n",
"\n",
"In this example, we perform optimization for a robust single qubit Hadamard gate of a qutrit system while [minimizing leakage out of the computational subspace](https://docs.q-ctrl.com/boulder-opal/user-guides/how-to-create-leakage-robust-single-qubit-gates).\n",
"The system is described by the following Hamiltonian:\n",
"$$\n",
"H(t) = \\frac{\\chi}{2} (a^\\dagger)^2 a^2 + (1+\\beta(t)) \\left(\\gamma(t) a + \\gamma^*(t) a^\\dagger \\right) ,\n",
"$$\n",
"where $\\chi$ is the anharmonicity, $\\gamma(t)$ is a complex time-dependent pulse, $\\beta(t)$ is a small, slowly-varying stochastic amplitude noise process, and $a = |0 \\rangle \\langle 1 | + \\sqrt{2} |1 \\rangle \\langle 2 |$.\n",
"\n",
"In this example, we will parametrize the real and imaginary parts of the pulse, $\\gamma(t) = \\gamma_I(t) + i \\gamma_Q(t)$, as a superposition of Lorentz functions,\n",
"$$\n",
"\\gamma_{I(Q)}(t) = \\sum_{n=1}^N c^{I(Q)}_n \\frac{\\sigma^2}{(t - t_n)^2 + \\sigma^2} ,\n",
"$$\n",
"where $c^{I(Q)}_n$ are the different real-valued coefficients describing the parametrization, $\\{t_n\\}$ are the centers of the Lorentz functions, and $\\sigma$ is their width.\n",
"This is a good choice for implementation in bandwidth-limited hardware as it is composed of smooth functions that go to zero away from the choice of centers $\\{t_n\\}$.\n",
"\n",
"Note that you can create a wide variety of analytical functions and superpositions using the different [mathematical and arithmetic graph operations](https://docs.q-ctrl.com/references/qctrl/Graphs.html#mathematical-functions)."
]
},
{
"cell_type": "code",
"execution_count": 1,
"metadata": {
"pycharm": {
"is_executing": true
}
},
"outputs": [],
"source": [
"import matplotlib.pyplot as plt\n",
"import numpy as np\n",
"\n",
"from qctrlvisualizer import get_qctrl_style, plot_controls\n",
"\n",
"plt.style.use(get_qctrl_style())\n",
"\n",
"from qctrl import Qctrl\n",
"\n",
"# Start a Boulder Opal session.\n",
"qctrl = Qctrl()"
]
},
{
"cell_type": "code",
"execution_count": 2,
"metadata": {
"pycharm": {
"is_executing": true
}
},
"outputs": [
{
"data": {
"application/vnd.jupyter.widget-view+json": {
"model_id": "",
"version_major": 2,
"version_minor": 0
},
"text/plain": [
" 0%| | 0/100 [00:00"
]
},
"metadata": {},
"output_type": "display_data"
}
],
"source": [
"# Define the function to generate the pulse components.\n",
"def custom_optimizable_superposition(graph, duration, coefficient_count, name):\n",
" \"\"\"\n",
" Create an STF optimizable superposition of Lorentzian functions.\n",
"\n",
" Parameters\n",
" ----------\n",
" graph : The graph where the signal will belong.\n",
" duration : The duration of the signal.\n",
" coefficient_count : The number of terms in the superposition.\n",
" name : The name of the Tensor node with the optimizable coefficients.\n",
"\n",
" Returns\n",
" -------\n",
" Stf\n",
" An optimizable superposition of Lorentzian functions.\n",
" \"\"\"\n",
"\n",
" # Define optimizable coefficients.\n",
" coefficients = graph.optimization_variable(\n",
" coefficient_count, lower_bound=-1, upper_bound=1, name=name\n",
" )\n",
"\n",
" # Define Lorentz function parameters.\n",
" width = 0.1 * duration\n",
" centers = np.linspace(0.3, 0.7, coefficient_count) * duration\n",
"\n",
" # Create Lorentz function superposition.\n",
" time = graph.identity_stf()\n",
" return graph.stf_sum(\n",
" [\n",
" coefficients[index] * width**2 / ((time - center) ** 2 + width**2)\n",
" for index, center in zip(range(coefficient_count), centers)\n",
" ]\n",
" )\n",
"\n",
"\n",
"# Define target and projector matrices.\n",
"hadamard = np.array([[1.0, 1.0, 0], [1.0, -1.0, 0], [0, 0, np.sqrt(2)]]) / np.sqrt(2)\n",
"qubit_projector = np.diag([1.0, 0.0, 0.0])\n",
"\n",
"# Define physical constraints\n",
"chi = -2 * np.pi * 300e6 # Hz\n",
"gamma_max = 2 * np.pi * 90e6 # Hz\n",
"segment_count = 200\n",
"duration = 100e-9 # s\n",
"sample_times = np.linspace(0, duration, segment_count)\n",
"coefficient_count = 5\n",
"\n",
"# Create graph object.\n",
"graph = qctrl.create_graph()\n",
"\n",
"# Define standard matrices.\n",
"a = graph.annihilation_operator(3)\n",
"ada = graph.number_operator(3)\n",
"\n",
"# Create gamma(t) signal using custom basis.\n",
"gamma_i = custom_optimizable_superposition(\n",
" graph, duration, coefficient_count, name=\"gamma_i\"\n",
")\n",
"gamma_q = custom_optimizable_superposition(\n",
" graph, duration, coefficient_count, name=\"gamma_q\"\n",
")\n",
"gamma = gamma_max * (gamma_i + 1j * gamma_q)\n",
"\n",
"# Discretize gamma to export and plot.\n",
"sample_gamma = graph.discretize_stf(\n",
" stf=gamma, duration=duration, segment_count=segment_count, name=r\"$\\gamma$\"\n",
")\n",
"\n",
"# Create anharmonicity term.\n",
"anharmonicity = chi / 2 * (ada @ ada - ada)\n",
"\n",
"# Create drive term.\n",
"drive = graph.hermitian_part(2 * gamma * a)\n",
"\n",
"# Create target operator in qubit subspace.\n",
"target_operator = graph.target(\n",
" hadamard.dot(qubit_projector), filter_function_projector=qubit_projector\n",
")\n",
"\n",
"# Create infidelity.\n",
"infidelity = graph.infidelity_stf(\n",
" hamiltonian=anharmonicity + drive,\n",
" target=target_operator,\n",
" sample_times=sample_times,\n",
" noise_operators=[drive],\n",
" name=\"infidelity\",\n",
")\n",
"\n",
"# Run the optimization.\n",
"optimization_result = qctrl.functions.calculate_optimization(\n",
" graph=graph,\n",
" cost_node_name=\"infidelity\",\n",
" output_node_names=[r\"$\\gamma$\", \"gamma_i\", \"gamma_q\"],\n",
" optimization_count=4,\n",
")\n",
"\n",
"# Retrieve results and plot optimized pulse.\n",
"coefficients_i = gamma_max * optimization_result.output.pop(\"gamma_q\")[\"value\"]\n",
"coefficients_q = gamma_max * optimization_result.output.pop(\"gamma_i\")[\"value\"]\n",
"\n",
"print(f\"\\nOptimized cost: {optimization_result.cost:.3e}\")\n",
"print(\"Optimized coefficients:\")\n",
"np.set_printoptions(precision=3)\n",
"print(\"\\tI:\", coefficients_i)\n",
"print(\"\\tQ:\", coefficients_q)\n",
"\n",
"plot_controls(optimization_result.output, smooth=True, polar=False)\n",
"plt.suptitle(\"Optimized control\")\n",
"plt.show()"
]
},
{
"cell_type": "code",
"execution_count": 3,
"metadata": {},
"outputs": [
{
"name": "stdout",
"output_type": "stream",
"text": [
"| Package | Version |\n",
"| --------------------- | ------------ |\n",
"| Python | 3.10.8 |\n",
"| matplotlib | 3.6.3 |\n",
"| numpy | 1.24.1 |\n",
"| scipy | 1.10.0 |\n",
"| qctrl | 20.1.1 |\n",
"| qctrl-commons | 17.7.0 |\n",
"| boulder-opal-toolkits | 2.0.0-beta.3 |\n",
"| qctrl-visualizer | 4.4.0 |\n"
]
}
],
"source": [
"from qctrl.utils import print_environment_related_packages\n",
"\n",
"print_environment_related_packages()"
]
}
],
"metadata": {
"kernelspec": {
"display_name": "Python 3 (ipykernel)",
"language": "python",
"name": "python3"
},
"language_info": {
"codemirror_mode": {
"name": "ipython",
"version": 3
},
"file_extension": ".py",
"mimetype": "text/x-python",
"name": "python",
"nbconvert_exporter": "python",
"pygments_lexer": "ipython3",
"version": "3.10.8"
}
},
"nbformat": 4,
"nbformat_minor": 4
}