diff --git a/doc/source/code-examples/applications-by-algorithm.rst b/doc/source/code-examples/applications-by-algorithm.rst index 3819ba99b3..1f69431d96 100644 --- a/doc/source/code-examples/applications-by-algorithm.rst +++ b/doc/source/code-examples/applications-by-algorithm.rst @@ -76,4 +76,12 @@ Diagonalization Algorithms .. toctree:: :maxdepth: 1 - tutorials/dbi/dbi.ipynb + tutorials/dbi/README.md + + tutorials/dbi/dbi_cost_functions.ipynb + tutorials/dbi/dbi_gradient_descent_strategies.ipynb + tutorials/dbi/dbi_scheduling.ipynb + tutorials/dbi/dbi_strategies_compare.ipynb + tutorials/dbi/dbi_strategy_Ising_model.ipynb + tutorials/dbi/dbi_strategy_Pauli-Z.ipynb + tutorials/dbi/dbi_tutorial_basic_intro.ipynb diff --git a/doc/source/code-examples/applications-by-topic.rst b/doc/source/code-examples/applications-by-topic.rst index 8d44557d48..500300e9b4 100644 --- a/doc/source/code-examples/applications-by-topic.rst +++ b/doc/source/code-examples/applications-by-topic.rst @@ -61,7 +61,6 @@ Quantum Physics tutorials/bell-variational/README.md tutorials/falqon/README.md tutorials/grover/README.md - tutorials/dbi/dbi.ipynb Quantum Machine Learning ^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/source/code-examples/applications.rst b/doc/source/code-examples/applications.rst index 8787659178..9fa6976b54 100644 --- a/doc/source/code-examples/applications.rst +++ b/doc/source/code-examples/applications.rst @@ -69,8 +69,6 @@ Quantum Physics tutorials/bell-variational/README.md tutorials/falqon/README.md tutorials/grover/README.md - tutorials/dbi/dbi.ipynb - Quantum Machine Learning ^^^^^^^^^^^^^^^^^^^^^^^^ diff --git a/doc/source/code-examples/tutorials/dbi/README.md b/doc/source/code-examples/tutorials/dbi/README.md index 78af756acd..50b9e9eaec 120000 --- a/doc/source/code-examples/tutorials/dbi/README.md +++ b/doc/source/code-examples/tutorials/dbi/README.md @@ -1,12 +1 @@ -# Double-bracket quantum algorithms - -Qibo features a model implementing double-bracke quantum algorithms (DBQAs) which are helpful for approximating eigenstates based on the ability to run the evolution under the input Hamiltonian. - -More specifically, given a Hamiltonian $H_0$, how can we find a circuit which after applying to the reference state (usually $|0\rangle^{\otimes L}$ for $L$ qubits) will approximate an eigenstate? - -A standard way is to run variational quantum circuits. For example, Qibo already features the `VQE` model [2] which provides the implementation of the variational quantum eigensolver framework. -DBQAs allow to go beyond VQE in that they take a different approach to compiling the quantum circuit approximating the eigenstate. - -[1] https://arxiv.org/abs/2206.11772 - -[2] https://github.com/qiboteam/vqe-sun +../../../../../examples/dbi/README.md \ No newline at end of file diff --git a/doc/source/code-examples/tutorials/dbi/dbi.ipynb b/doc/source/code-examples/tutorials/dbi/dbi.ipynb deleted file mode 120000 index 7deb426d68..0000000000 --- a/doc/source/code-examples/tutorials/dbi/dbi.ipynb +++ /dev/null @@ -1 +0,0 @@ -../../../../../examples/dbi/dbi.ipynb \ No newline at end of file diff --git a/doc/source/code-examples/tutorials/dbi/dbi_cost_functions.ipynb b/doc/source/code-examples/tutorials/dbi/dbi_cost_functions.ipynb new file mode 120000 index 0000000000..1c36a4d760 --- /dev/null +++ b/doc/source/code-examples/tutorials/dbi/dbi_cost_functions.ipynb @@ -0,0 +1 @@ +../../../../../examples/dbi/dbi_cost_functions.ipynb \ No newline at end of file diff --git a/doc/source/code-examples/tutorials/dbi/dbi_strategies_compare.ipynb b/doc/source/code-examples/tutorials/dbi/dbi_strategies_compare.ipynb new file mode 120000 index 0000000000..78936e1a6b --- /dev/null +++ b/doc/source/code-examples/tutorials/dbi/dbi_strategies_compare.ipynb @@ -0,0 +1 @@ +../../../../../examples/dbi/dbi_strategies_compare.ipynb \ No newline at end of file diff --git a/examples/dbi/README.md b/examples/dbi/README.md new file mode 100644 index 0000000000..76640ff99f --- /dev/null +++ b/examples/dbi/README.md @@ -0,0 +1,58 @@ +# Double-bracket quantum algorithms + +Qibo features a model implementing double-bracke quantum algorithms (DBQAs) which are helpful for approximating eigenstates based on the ability to run the evolution under the input Hamiltonian. + +More specifically, given a Hamiltonian $H_0$, how can we find a circuit which after applying to the reference state (usually $|0\rangle^{\otimes L}$ for $L$ qubits) will approximate an eigenstate? + +A standard way is to run variational quantum circuits. For example, Qibo already features the `VQE` model [2] which provides the implementation of the variational quantum eigensolver framework. +DBQAs allow to go beyond VQE in that they take a different approach to compiling the quantum circuit approximating the eigenstate. + +## What is the unitary of DBQA? + +Given $H_0$ we begin by assuming that we were given a diagonal and hermitian operator $D_0$ and a time $s_0$. +The `dbi` module provides numerical strategies for selecting them. +For any such choice we define the bracket +$$ W_0 = [D_0, H_0]$$ +and the double-bracket rotation (DBR) of the input Hamiltonian to time $s$ +$$H_0(s) = e^{sW} H e^{- s W}$$ + +### Why are double-bracket rotations useful? +We can show that the magnitude of the off-diagonal norms will decrease. +For this let us set the notation that $\sigma(A)$ is the restriction to the off-diagonal of the matrix A. +In `numpy` this can be implemented by `\sigma(A) = A-np.diag(A)`. In Qibo we implement this as +https://github.com/qiboteam/qibo/blob/8c9c610f5f2190b243dc9120a518a7612709bdbc/src/qibo/models/dbi/double_bracket.py#L145-L147 +which is part of the basic `DoubleBracketIteration` class in the `dbi` module. + +With this notation we next use the Hilbert-Schmidt scalar product and norm to measure the progress of diagonalization + $$||\sigma(H_0(s))||^2- ||\sigma (H_0 )||^2= -2s \langle W, [H,\sigma(H)\rangle+O(s^2)$$ +This equation tells us that as long as the scalar product $\langle W, [H,\sigma(H)\rangle$ is positive then after the DBR the magnitude of the off-diagonal couplings in $H_0(s)$ is less than in $H_0$. + +For the implementation of the DBR unitary $U_0(s) = e^{-s W_0}$ see +https://github.com/qiboteam/qibo/blob/363a6e5e689e5b907a7602bd1cc8d9811c60ee69/src/qibo/models/dbi/double_bracket.py#L68 + +### How to choose $D$? + +For theoretical considerations the canonical bracket is useful. +For this we need the notation of the dephasing channel $\Delta(H)$ which is equivalent to `np.diag(h)`. + $M = [\Delta(H),\sigma(H)]= [H,\sigma(H)]= [\Delta(H),H]$ + The canonical bracket appears on its own in the monotonicity relation above and gives an unconditional reduction of the magnitude of the off-diagonal terms + $$||\sigma(H_0(s))||^2- ||\sigma (H_0 )||^2= -2s ||M||^2+O(s^2)$$ +- the multi qubit Pauli Z generator with $Z(\mu) = (Z_1)^{\mu_1}\ldots (Z_L)^{\mu_L}$ where we optimize over all binary strings $\mu\in \{0,1\}^L$ +- the magnetic field $D = \sum_i B_i Z_i$ +- the two qubit Ising model $D = \sum_i B_i Z_i + \sum_{i,j} J_{i,j} Z_i Z_j$, please follow the tutorial by Matteo and use the QIBO ising model for that with $h=0$ + + +### How to choose s? + +The theory above shows that in generic cases the DBR will have a linear diagonalization effect (as quantified by $||\sigma(H_0(s))||$). +This can be further expanded with Taylor expansion and the Qibo implementation comes with methods for fitting the first local minimum. +Additionally a grid search for the optimal step is provided for an exhaustive evaluation and hyperopt can be used for a more efficient 'unstructured' optimization; additionally simulated annealing is provided which sometimes outperforms hyperopt (and grid search), see example notebooks. +The latter methods may output DBR durations $s_k$ which correspond to secondary local minima. + + + + + +[1] https://arxiv.org/abs/2206.11772 + +[2] https://github.com/qiboteam/vqe-sun diff --git a/examples/dbi/dbi_cost_functions.ipynb b/examples/dbi/dbi_cost_functions.ipynb new file mode 100644 index 0000000000..015a441a94 --- /dev/null +++ b/examples/dbi/dbi_cost_functions.ipynb @@ -0,0 +1,349 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Double-bracket Iteration other cost functions\n", + "\n", + "This notebook presents two additional cost functions for the double-bracket flow: least-squares and energy fluctuation with their respectice scheduling methods." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from copy import deepcopy\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from qibo import hamiltonians, set_backend\n", + "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketScheduling, DoubleBracketIteration, DoubleBracketCostFunction" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Least-squares\n", + "\n", + "The cost function is defined as: $\\frac{1}{2}||D-H_k||^2 =\\frac{1}{2}(||D||^2+||H||^2) -Tr(D H_k)$ as in (the negative of https://epubs.siam.org/doi/abs/10.1137/S0036141092229732?journalCode=sjmael) We seek to minimize this function at each DBF iteration. For numerical optimizations, we also ignore the norm of H term as for a given hamiltonian it is fixed through out the flow.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Hamiltonian\n", + "set_backend(\"numpy\")\n", + "\n", + "# hamiltonian parameters\n", + "nqubits = 5\n", + "h = 3.0\n", + "\n", + "# define the hamiltonian\n", + "H_TFIM = hamiltonians.TFIM(nqubits=nqubits, h=h)\n", + "\n", + "# define the least-squares cost function\n", + "cost = DoubleBracketCostFunction.least_squares\n", + "# initialize class\n", + "dbi = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.single_commutator,cost=cost)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# generate data for plotting sigma decrease of the first step\n", + "d = np.diag(np.linspace(1,2**nqubits,2**nqubits))/2**nqubits\n", + "s_space = np.linspace(1e-5, 1.0, 500)\n", + "off_diagonal_norm_diff = []\n", + "potential = []\n", + "for s in s_space:\n", + " dbi_eval = deepcopy(dbi)\n", + " dbi_eval(s,d=d)\n", + " off_diagonal_norm_diff.append(dbi_eval.off_diagonal_norm - dbi.off_diagonal_norm)\n", + " potential.append(dbi_eval.least_squares(d=d))\n", + "\n", + "# grid_search\n", + "step_grid = dbi.choose_step(scheduling=DoubleBracketScheduling.grid_search,d=d)\n", + "print('grid_search step:', step_grid)\n", + "# hyperopt\n", + "step_hyperopt = dbi.choose_step(scheduling=DoubleBracketScheduling.hyperopt,d=d, max_evals=100, step_max=0.6)\n", + "print('hyperopt_search step:', step_hyperopt)\n", + "# polynomial\n", + "step_poly = dbi.choose_step(scheduling=DoubleBracketScheduling.polynomial_approximation,d=d, n=3)\n", + "print('polynomial_approximation step:', step_poly)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the results\n", + "plt.figure()\n", + "plt.plot(s_space, potential)\n", + "plt.xlabel('s')\n", + "plt.axvline(x=step_grid, color='r', linestyle='-',label='grid_search')\n", + "plt.axvline(x=step_hyperopt, color='g', linestyle='--',label='hyperopt')\n", + "plt.axvline(x=step_poly, color='m', linestyle='-.',label='polynomial')\n", + "plt.title('First DBI step')\n", + "plt.ylabel('Least squares cost function')\n", + "plt.legend()\n", + "plt.figure()\n", + "plt.plot(s_space, off_diagonal_norm_diff)\n", + "plt.axvline(x=step_grid, color='r', linestyle='-',label='grid_search')\n", + "plt.axvline(x=step_hyperopt, color='g', linestyle='--',label='hyperopt')\n", + "plt.axvline(x=step_poly, color='m', linestyle='-.',label='polynomial')\n", + "plt.ylabel(r'$||\\sigma(H_0)||-\\sigma(H_k)||$')\n", + "plt.xlabel('s')\n", + "plt.title('First DBI step')\n", + "plt.legend()\n", + "print('The minimum for cost function in the tested range is:', step_grid)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Comparison of the least-squares cost function with the original cost function using the polynomial scheduling method" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "d = np.diag(np.linspace(1,2**nqubits,2**nqubits))\n", + "off_diagonal_norm_diff = [dbi.off_diagonal_norm]\n", + "off_diagonal_norm_diff_least_squares = [dbi.off_diagonal_norm]\n", + "iters = 100\n", + "dbi_ls = deepcopy(dbi)\n", + "cost = DoubleBracketCostFunction.off_diagonal_norm\n", + "dbi_od = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.single_commutator,cost=cost)\n", + "for _ in range(iters):\n", + " step_poly = dbi_od.choose_step(scheduling=DoubleBracketScheduling.polynomial_approximation, d=d, n=3)\n", + " dbi_od(step_poly,d=d)\n", + " step_poly = dbi_ls.choose_step(scheduling=DoubleBracketScheduling.polynomial_approximation, d=d, n=3)\n", + " dbi_ls(step_poly,d=d)\n", + " off_diagonal_norm_diff.append(dbi_od.off_diagonal_norm)\n", + " off_diagonal_norm_diff_least_squares.append(dbi_ls.off_diagonal_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.plot(range(iters+1), off_diagonal_norm_diff, label=r'Off-diagonal norm')\n", + "plt.plot(range(iters+1), off_diagonal_norm_diff_least_squares, label=r'Least squares')\n", + "plt.xlabel('Iterations')\n", + "plt.ylabel(r'$||\\sigma(H_k)||$')\n", + "plt.legend()\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Energy fluctuation\n", + "\n", + "This cost function is defined as: $\\Xi_k^2 (\\mu) = \\langle \\mu | H_k^2| \\mu \\rangle - \\langle \\mu | H_k| \\mu \\rangle^2$. We must specify the state $| \\mu \\rangle$ for which we want to minimize the fluctuation. The overall diagonalization isn't guaranteed.\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Hamiltonian\n", + "set_backend(\"numpy\")\n", + "\n", + "# hamiltonian parameters\n", + "nqubits = 3\n", + "h = 3.0\n", + "\n", + "# define the hamiltonian\n", + "H_TFIM = hamiltonians.TFIM(nqubits=nqubits, h=h)\n", + "\n", + "# define the energy fluctuation cost function\n", + "cost = DoubleBracketCostFunction.off_diagonal_norm\n", + "# define the state\n", + "state = np.zeros(2**nqubits)\n", + "state[3] = 1\n", + "# initialize class\n", + "dbi = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.single_commutator,cost=cost, ref_state=state)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# generate data for plotting sigma decrease of the first step\n", + "d = np.diag(np.linspace(2**nqubits,1,2**nqubits))/2**nqubits\n", + "s_space = np.linspace(-1, 1, 1000)\n", + "off_diagonal_norm_diff = []\n", + "fluctuation = []\n", + "for s in s_space:\n", + " dbi_eval = deepcopy(dbi)\n", + " dbi_eval(s,d=d)\n", + " off_diagonal_norm_diff.append(dbi_eval.off_diagonal_norm - dbi.off_diagonal_norm)\n", + " fluctuation.append(dbi_eval.energy_fluctuation(state=state))\n", + "\n", + "# grid_search\n", + "step_grid = dbi.choose_step(scheduling=DoubleBracketScheduling.grid_search,d=d)\n", + "print('grid_search step:', step_grid)\n", + "# hyperopt\n", + "step_hyperopt = dbi.choose_step(scheduling=DoubleBracketScheduling.hyperopt,d=d, max_evals=100, step_max=0.6)\n", + "print('hyperopt_search step:', step_hyperopt)\n", + "# polynomial\n", + "step_poly = dbi.choose_step(scheduling=DoubleBracketScheduling.polynomial_approximation,d=d, n=3)\n", + "print('polynomial_approximation step:', step_poly)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Plot the results\n", + "plt.figure()\n", + "plt.plot(s_space, fluctuation)\n", + "plt.xlabel('s')\n", + "plt.axvline(x=step_grid, color='r', linestyle='-',label='grid_search')\n", + "plt.axvline(x=step_hyperopt, color='g', linestyle='--',label ='hyperopt')\n", + "plt.axvline(x=step_poly, color='m', linestyle='-.',label='polynomial')\n", + "plt.title('First DBI step')\n", + "plt.ylabel('Energy fluctuation')\n", + "plt.legend()\n", + "plt.figure()\n", + "plt.plot(s_space, off_diagonal_norm_diff)\n", + "plt.axvline(x=step_grid, color='r', linestyle='-',label='grid_search')\n", + "plt.axvline(x=step_hyperopt, color='g', linestyle='--',label='hyperopt')\n", + "plt.axvline(x=step_poly, color='m', linestyle='-.',label='polynomial')\n", + "plt.ylabel(r'$||\\sigma(H_0)||-\\sigma(H_k)||$')\n", + "plt.xlabel('s')\n", + "plt.title('First DBI step')\n", + "plt.legend()\n", + "print('The minimum for cost function in the tested range is:', step_grid)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "d = np.diag(np.linspace(1,2**nqubits,2**nqubits))\n", + "off_diagonal_norm_diff = [dbi.off_diagonal_norm]\n", + "energy_fluc = [dbi.energy_fluctuation(state=state)]\n", + "iters = 10\n", + "dbi_ = deepcopy(dbi)\n", + "for _ in range(iters):\n", + " step_poly = dbi_.choose_step(scheduling=DoubleBracketScheduling.polynomial_approximation, d=d, n=3)\n", + " dbi_(step_poly,d=d)\n", + " off_diagonal_norm_diff.append(dbi_.off_diagonal_norm)\n", + " energy_fluc.append(dbi_.energy_fluctuation(state=state))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "plt.plot(range(iters+1), off_diagonal_norm_diff)\n", + "plt.xlabel('Iterations')\n", + "plt.ylabel(r'$||\\sigma(H_k)||$')\n", + "\n", + "plt.figure()\n", + "plt.plot(range(iters+1), energy_fluc)\n", + "plt.xlabel('Iterations')\n", + "plt.ylabel(r'Energy fluctuation')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iters = 30\n", + "states = [0,1,2,3,4,5,6,7]\n", + "energy = np.empty((len(states),iters))\n", + "\n", + "\n", + "d = (np.diag(np.linspace(1,2**nqubits,2**nqubits)))\n", + "for i in range(len(states)):\n", + " dbi_ = deepcopy(dbi)\n", + " dbi_.state = states[i]\n", + " for j in range(iters):\n", + " step_poly = dbi_.choose_step(scheduling=DoubleBracketScheduling.polynomial_approximation, d=d, n=3)\n", + " if step_poly is not None:\n", + " dbi_(step_poly, d=d)\n", + " energy[i,j] = np.real(dbi_.h.matrix[states[i],states[i]])\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "eigvals = np.linalg.eigh(dbi_.h.matrix)[0]\n", + "print('Eigenvalues:', eigvals )\n", + "plt.figure()\n", + "for i in range(len(states)):\n", + " plt.plot(range(iters), energy[i,:],'.', label='State ' + str(states[i]))\n", + "for eigvals in eigvals:\n", + " plt.axhline(y=eigvals, color='r', linestyle='--')\n", + "plt.xlabel('Iterations')\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "test_qibo", + "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.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/dbi/dbi_gradient_descent_strategies.ipynb b/examples/dbi/dbi_gradient_descent_strategies.ipynb new file mode 100644 index 0000000000..58b6ec470b --- /dev/null +++ b/examples/dbi/dbi_gradient_descent_strategies.ipynb @@ -0,0 +1,343 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Double-bracket Gradient Descent Stratgies\n", + "This notebook demonstrates the gradient descent strategies for double-bracket rotations. The mehods uses a numerical method to find the gradient of the cost function with respect to the diagonal operator, and thereby variate the diagonal operator of the rotation. \n", + "\n", + "Finding the gradient requires the parameterization of the diagonal operator, and there are two ways of doing so:\n", + "\n", + "1. Pauli-basis: $D(B,J)= \\sum B_i Z_i + \\sum J_{ij}Z_iZ_j + ...$\n", + "2. Computational-basis: $D(A)=\\sum A_i|i\\rangle\\langle i|$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from qibo.models.dbi.double_bracket import DoubleBracketIteration, DoubleBracketGeneratorType, DoubleBracketScheduling, DoubleBracketCostFunction\n", + "from qibo.models.dbi.utils import generate_pauli_operator_dict, decompose_into_pauli_basis, params_to_diagonal_operator, ParameterizationTypes\n", + "from copy import deepcopy\n", + "from qibo.models.dbi.utils_dbr_strategies import gradient_descent\n", + "import numpy as np\n", + "from qibo import set_backend, hamiltonians\n", + "from qibo.hamiltonians import Hamiltonian\n", + "from qibo.quantum_info import random_hermitian\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def visualize_matrix(matrix, title=\"\"):\n", + " \"\"\"Visualize hamiltonian in a heatmap form.\"\"\"\n", + " fig, ax = plt.subplots(figsize=(5,5))\n", + " ax.set_title(title)\n", + " try:\n", + " im = ax.imshow(np.absolute(matrix), cmap=\"inferno\")\n", + " except TypeError:\n", + " im = ax.imshow(np.absolute(matrix.get()), cmap=\"inferno\")\n", + " fig.colorbar(im, ax=ax)\n", + "\n", + "def s_hist_to_plot(s_hist):\n", + " # convert list of step durations taken to plotable\n", + " s_plot = [0] * len(s_hist)\n", + " for i in range(len(s_hist)):\n", + " if i != 0:\n", + " s_plot[i] = s_plot[i-1] + s_hist[i]\n", + " return s_plot" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Random Hamiltonian" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# set the qibo backend (we suggest qibojit if N >= 20)\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", + "\n", + "# hamiltonian parameters\n", + "nqubits = 5\n", + "seed = 10\n", + "\n", + "# define the hamiltonian\n", + "h0 = random_hermitian(2**nqubits, seed=seed)\n", + "dbi = DoubleBracketIteration(\n", + " Hamiltonian(nqubits, h0),\n", + " mode=DoubleBracketGeneratorType.single_commutator,\n", + " scheduling=DoubleBracketScheduling.hyperopt,\n", + " cost=DoubleBracketCostFunction.off_diagonal_norm\n", + ")\n", + "# vosualize the matrix\n", + "visualize_matrix(dbi.h.matrix, title=\"Target hamiltonian\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Then we set up the required parameters for gradient descent." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pauli-basis\n", + "pauli_operator_dict = generate_pauli_operator_dict(nqubits)\n", + "pauli_operators = list(pauli_operator_dict.values())\n", + "# let initial d be approximation of $\\Delta(H)\n", + "d_coef_pauli = decompose_into_pauli_basis(dbi.diagonal_h_matrix, pauli_operators=pauli_operators)\n", + "d_pauli = sum([d_coef_pauli[i]*pauli_operators[i] for i in range(nqubits)])\n", + "\n", + "# Computational basis\n", + "d_coef_computational_partial = d_pauli.diagonal()\n", + "d_coef_computational_full = dbi.diagonal_h_matrix.diagonal()\n", + "d_computational_partial = params_to_diagonal_operator(d_coef_computational_partial, nqubits, ParameterizationTypes.computational, normalize=False)\n", + "d_computational_full = params_to_diagonal_operator(d_coef_computational_full, nqubits, ParameterizationTypes.computational, normalize=False)\n", + "\n", + "plt.plot(d_coef_computational_partial, label=\"computational basis partial\")\n", + "plt.plot(d_coef_computational_full, label=r\"computational basis full = $\\Delta(H)$\")\n", + "plt.legend()\n", + "plt.title(r\"Diagonal entries of $D$\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Now we want to compare 3 scenarios:\n", + "\n", + "1. Pauli-basis: an approximation to the diagonal of $H$\n", + "2. Computational-partial: same as 1. in the computational basis.\n", + "3. Computational-full: a full parameterization of the diagonal of $H$ in the computational basis." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 1. Pauli-basis\n", + "NSTEPS = 5\n", + "dbi_pauli = deepcopy(dbi)\n", + "loss_hist_pauli, d_params_hist_pauli, s_hist_pauli = gradient_descent(dbi_pauli, NSTEPS, d_coef_pauli, ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 2. Computational_partial\n", + "dbi_computational_partial = deepcopy(dbi)\n", + "loss_hist_computational_partial, d_params_hist_computational_partiali, s_computational_partial = gradient_descent(dbi_computational_partial, NSTEPS, d_coef_computational_partial, ParameterizationTypes.computational)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 3. Computational_full\n", + "dbi_computational_full = deepcopy(dbi)\n", + "loss_hist_computational_full, d_params_hist_computational_full, s_computational_full = gradient_descent(dbi_computational_full, NSTEPS, d_coef_computational_full, ParameterizationTypes.computational)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s_plot_pauli = s_hist_to_plot(s_hist_pauli)\n", + "s_plot_computational_partial = s_hist_to_plot(s_computational_partial)\n", + "s_plot_computational_full = s_hist_to_plot(s_computational_full)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(s_plot_pauli, loss_hist_pauli, label=\"pauli basis\", marker=\"o\")\n", + "plt.plot(s_plot_computational_partial, loss_hist_computational_partial, label=\"computational partial\", marker=\"o\")\n", + "plt.plot(s_plot_computational_full, loss_hist_computational_full, label=\"computational full\", marker=\"o\")\n", + "plt.legend()\n", + "plt.title(\"Off-diagonal norm\")\n", + "plt.ylabel(r\"$||\\sigma(H)||_{HS}$\")\n", + "plt.xlabel(\"s\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# TFIM" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# hamiltonian parameters\n", + "nqubits = 5\n", + "h = 3\n", + "\n", + "# define the hamiltonian\n", + "h = hamiltonians.TFIM(nqubits=nqubits, h=h)\n", + "dbi = DoubleBracketIteration(\n", + " h,\n", + " mode=DoubleBracketGeneratorType.single_commutator,\n", + " scheduling=DoubleBracketScheduling.hyperopt\n", + ")\n", + "# vosualize the matrix\n", + "visualize_matrix(dbi.h.matrix, title=\"Target hamiltonian\")" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Pauli-basis\n", + "pauli_operator_dict = generate_pauli_operator_dict(nqubits)\n", + "pauli_operators = list(pauli_operator_dict.values())\n", + "# let initial d be approximation of $\\Delta(H)\n", + "d_coef_pauli = decompose_into_pauli_basis(dbi.diagonal_h_matrix, pauli_operators=pauli_operators)\n", + "d_pauli = sum([d_coef_pauli[i]*pauli_operators[i] for i in range(nqubits)])\n", + "\n", + "# Computational basis\n", + "d_coef_computational_partial = d_pauli.diagonal()\n", + "d_coef_computational_full = dbi.diagonal_h_matrix.diagonal()\n", + "d_computational_partial = params_to_diagonal_operator(d_coef_computational_partial, nqubits, ParameterizationTypes.computational, normalize=False)\n", + "d_computational_full = params_to_diagonal_operator(d_coef_computational_full, nqubits, ParameterizationTypes.computational, normalize=False)\n", + "\n", + "plt.plot(d_coef_computational_partial, label=\"computational basis partial\")\n", + "plt.plot(d_coef_computational_full, label=r\"computational basis full = $\\Delta(H)$\")\n", + "plt.legend()\n", + "plt.title(r\"Diagonal entries of $D$\")\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 1. Pauli-basis\n", + "NSTEPS = 3\n", + "dbi_pauli = deepcopy(dbi)\n", + "loss_hist_pauli, d_params_hist_pauli, s_hist_pauli = gradient_descent(dbi_pauli, NSTEPS, d_coef_pauli, ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 2. Computational_partial\n", + "dbi_computational_partial = deepcopy(dbi)\n", + "loss_hist_computational_partial, d_params_hist_computational_partiali, s_computational_partial = gradient_descent(dbi_computational_partial, NSTEPS, d_coef_computational_partial, ParameterizationTypes.computational)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# 3. Computational_full\n", + "dbi_computational_full = deepcopy(dbi)\n", + "loss_hist_computational_full, d_params_hist_computational_full, s_computational_full = gradient_descent(dbi_computational_full, NSTEPS, d_coef_computational_full, ParameterizationTypes.computational)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "s_plot_pauli = s_hist_to_plot(s_hist_pauli)\n", + "s_plot_computational_partial = s_hist_to_plot(s_computational_partial)\n", + "s_plot_computational_full = s_hist_to_plot(s_computational_full)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.plot(s_plot_pauli, loss_hist_pauli, label=\"pauli basis\", marker=\"o\")\n", + "plt.plot(s_plot_computational_partial, loss_hist_computational_partial, label=\"computational partial\", marker=\"o\")\n", + "plt.plot(s_plot_computational_full, loss_hist_computational_full, label=\"computational full\", marker=\"o\")\n", + "plt.legend()\n", + "plt.title(\"Off-diagonal norm\")\n", + "plt.ylabel(r\"$||\\sigma(H)||_{HS}$\")\n", + "plt.xlabel(\"s\")\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "After changing the cost function and scheduling method, we notice that quite consistently, the Pauli-based parameterization diagonalizes the hamiltonian the best, and for the first few iterations, the Computational-based partial (same initial operator as Pauli) performs very similarly, and diverges later on." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "nqubits = 3\n", + "pauli_operator_dict = generate_pauli_operator_dict(\n", + " nqubits, parameterization_order=1\n", + ")\n", + "params = [1, 2, 3]\n", + "operator_pauli = sum([\n", + " params[i] * list(pauli_operator_dict.values())[i] for i in range(nqubits)\n", + "])\n", + "assert (\n", + " operator_pauli\n", + " == params_to_diagonal_operator(\n", + " params, nqubits=nqubits, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict\n", + " )\n", + ").all()\n", + "operator_element = params_to_diagonal_operator(\n", + " params, nqubits=nqubits, parameterization=ParameterizationTypes.computational\n", + ")\n", + "assert (operator_element.diagonal() == params).all()" + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/dbi/dbi_scheduling.ipynb b/examples/dbi/dbi_scheduling.ipynb index 066d56640d..fb9d349ebf 100644 --- a/examples/dbi/dbi_scheduling.ipynb +++ b/examples/dbi/dbi_scheduling.ipynb @@ -29,7 +29,10 @@ "\n", "from qibo import hamiltonians, set_backend\n", "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketScheduling, DoubleBracketIteration\n", - "from qibo.models.dbi.utils import *" + "from qibo.hamiltonians import SymbolicHamiltonian\n", + "from qibo.models.dbi.utils import str_to_symbolic, generate_Z_operators\n", + "from qibo.models.dbi.utils_scheduling import polynomial_step\n", + "from qibo.models.dbi.utils_dbr_strategies import select_best_dbr_generator" ] }, { @@ -47,7 +50,7 @@ "outputs": [], "source": [ "# Hamiltonian\n", - "set_backend(\"qibojit\", \"numba\")\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", "\n", "# hamiltonian parameters\n", "nqubits = 5\n", @@ -65,7 +68,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "We first generate the relationship between the step duration and the off-diagoanl norm (loss function) for the first step of the iteration." + "We first run a sweep of step duration to map the off-diagonal norm in this range." ] }, { @@ -87,7 +90,8 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "The default scheduling strategy is grid search: `DoubleBracketScheduling.use_grid_serach`. This strategy specifies a list of step durations to test one by one and finds the one that maximizes the cost function (off-digonal norm of Hamiltonian)" + "The default scheduling strategy is grid search: `DoubleBracketScheduling.\n", + "grid_serach`. This strategy specifies a list of step durations to test one by one and finds the one that maximizes the cost function (off-digonal norm of Hamiltonian)" ] }, { @@ -97,14 +101,15 @@ "outputs": [], "source": [ "# grid_search\n", - "step_grid = dbi.choose_step(scheduling=DoubleBracketScheduling.use_grid_search)\n", + "step_grid = dbi.choose_step(scheduling=DoubleBracketScheduling.grid_search)\n", "print('grid_search step:', step_grid)\n", "# hyperopt\n", - "step_hyperopt = dbi.choose_step(scheduling=DoubleBracketScheduling.use_hyperopt, max_evals=100, step_max=0.6)\n", + "step_hyperopt = dbi.choose_step(scheduling=DoubleBracketScheduling.hyperopt, max_evals=100, step_max=0.6)\n", "print('hyperopt_search step:', step_hyperopt)\n", - "# polynomial expansion\n", - "step_poly = dbi.choose_step(scheduling=DoubleBracketScheduling.use_polynomial_approximation, n=5)\n", - "print('polynomial_approximation step:', step_poly)" + "step_poly = dbi.choose_step(scheduling=DoubleBracketScheduling.polynomial_approximation, n=5)\n", + "print('polynomial_approximation step:', step_poly)\n", + "step_sa = dbi.choose_step(scheduling=DoubleBracketScheduling.simulated_annealing)\n", + "print('simulated_annealing step:', step_sa)" ] }, { @@ -118,9 +123,10 @@ "plt.axvline(x=step_grid, color='r', linestyle='-',label='grid_search')\n", "plt.axvline(x=step_hyperopt, color='g', linestyle='--',label='hyperopt')\n", "plt.axvline(x=step_poly, color='m', linestyle='-.',label='polynomial')\n", + "plt.axvline(x=step_sa, color='b', linestyle=':',label='simulated annealing')\n", "plt.ylabel(r'$||\\sigma(H_0)||-\\sigma(H_k)||$')\n", "plt.xlabel('s')\n", - "plt.title('hyperopt first step')\n", + "plt.title('First DBI step')\n", "plt.legend()\n", "print('The minimum for cost function in the tested range is:', step_grid)" ] @@ -141,8 +147,11 @@ "outputs": [], "source": [ "# Generate the digaonal operators\n", - "Z_op = SymbolicHamiltonian(str_to_symbolic(\"Z\"*nqubits)).dense.matrix\n", - "ZI_op = SymbolicHamiltonian(str_to_symbolic(\"Z\"*(nqubits-1)+\"I\")).dense.matrix" + "Z_str = \"Z\"*nqubits\n", + "ZI_str = \"Z\"*(nqubits-1)+\"I\"\n", + "Z_op = SymbolicHamiltonian(str_to_symbolic(Z_str)).dense.matrix\n", + "ZI_op = SymbolicHamiltonian(str_to_symbolic(ZI_str)).dense.matrix\n", + "op_dict = {Z_str:Z_op, ZI_str: ZI_op}" ] }, { @@ -152,7 +161,8 @@ "outputs": [], "source": [ "dbi = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.single_commutator)\n", - "d = Z_op\n", + "d_str = ZI_str\n", + "d = op_dict[d_str]\n", "# generate data for plotting sigma decrease of the first step\n", "s_space = np.linspace(1e-5, 0.6, 100)\n", "off_diagonal_norm_diff = []\n", @@ -169,14 +179,21 @@ "outputs": [], "source": [ "# grid_search\n", - "step_grid = dbi.choose_step(scheduling=DoubleBracketScheduling.use_grid_search, step_max=0.6, d=d)\n", - "print('grid_search step:', step_grid)\n", + "step_grid = dbi.choose_step(scheduling=DoubleBracketScheduling.grid_search, step_max=0.6, d=d)\n", + "grid_min = dbi.loss(step=step_grid, d=d)-dbi.off_diagonal_norm\n", + "print('grid_search step:', step_grid, 'loss', grid_min)\n", "# hyperopt\n", - "step_hyperopt = dbi.choose_step(scheduling=DoubleBracketScheduling.use_hyperopt, d=d, max_evals=100, step_max=0.6)\n", - "print('hyperopt_search step:', step_hyperopt)\n", + "step_hyperopt = dbi.choose_step(scheduling=DoubleBracketScheduling.hyperopt, d=d, max_evals=100, step_max=0.6)\n", + "hyperopt_min = dbi.loss(step=step_hyperopt, d=d)-dbi.off_diagonal_norm\n", + "print('hyperopt_search step:', step_hyperopt, 'loss', hyperopt_min)\n", "# polynomial expansion\n", - "step_poly = dbi.choose_step(scheduling=DoubleBracketScheduling.use_polynomial_approximation, d=d, n=5)\n", - "print('polynomial_approximation step:', step_poly)" + "step_poly = dbi.choose_step(scheduling=DoubleBracketScheduling.polynomial_approximation, d=d, n=5)\n", + "poly_min = dbi.loss(step=step_poly, d=d)-dbi.off_diagonal_norm\n", + "print('polynomial_approximation step:', step_poly, 'loss', poly_min)\n", + "# simulated annealing\n", + "step_sa = dbi.choose_step(scheduling=DoubleBracketScheduling.simulated_annealing, d=d)\n", + "sa_min = dbi.loss(step=step_sa, d=d)-dbi.off_diagonal_norm\n", + "print('simulated_annealing step:', step_sa, 'loss', sa_min)" ] }, { @@ -188,11 +205,15 @@ "# Plot the results\n", "plt.plot(s_space, off_diagonal_norm_diff)\n", "plt.axvline(x=step_grid, color='r', linestyle='-',label='grid_search')\n", + "plt.text(x=step_grid, y=grid_min, s=f'grid min \\n{round(grid_min,3)}')\n", + "plt.text(x=step_poly, y=poly_min, s=f'poly min \\n{round(poly_min,3)}')\n", + "plt.text(x=step_sa, y=sa_min, s=f'sa min \\n{round(sa_min,3)}')\n", "plt.axvline(x=step_hyperopt, color='g', linestyle='--',label='hyperopt')\n", "plt.axvline(x=step_poly, color='m', linestyle='-.',label='polynomial')\n", + "plt.axvline(x=step_sa, color='b', linestyle=':',label='simulated annealing')\n", "plt.ylabel(r'$||\\sigma(H_0)||-\\sigma(H_k)||$')\n", "plt.xlabel('s')\n", - "plt.title('hyperopt first step')\n", + "plt.title(f'First DBI step with D={d_str}')\n", "plt.legend()\n", "print('The minimum for cost function in the tested range is:', step_grid)" ] @@ -208,7 +229,7 @@ "cell_type": "markdown", "metadata": {}, "source": [ - "## Use polynomial expansion as an restriction of hyperopt/grid range" + "## Use polynomial expansion as an restriction for hyperopt/grid range" ] }, { @@ -225,10 +246,10 @@ " step_min = step_poly - search_range/2\n", " step_max = step_poly + search_range/2\n", "# grid_search\n", - "step_grid = dbi.choose_step(scheduling=DoubleBracketScheduling.use_grid_search, step_min=step_min, step_max=step_max, d=d)\n", + "step_grid = dbi.choose_step(scheduling=DoubleBracketScheduling.grid_search, step_min=step_min, step_max=step_max, d=d)\n", "print('grid_search step:', step_grid)\n", "# hyperopt\n", - "step_hyperopt = dbi.choose_step(scheduling=DoubleBracketScheduling.use_hyperopt, step_min=step_min, step_max=step_max, max_evals=100, d=d,)\n", + "step_hyperopt = dbi.choose_step(scheduling=DoubleBracketScheduling.hyperopt, step_min=step_min, step_max=step_max, max_evals=100, d=d,)\n", "print('hyperopt_search step:', step_hyperopt)" ] }, @@ -245,31 +266,190 @@ "plt.axvline(x=step_poly, color='m', linestyle='-.',label='polynomial')\n", "plt.ylabel(r'$||\\sigma(H_0)||-\\sigma(H_k)||$')\n", "plt.xlabel('s')\n", - "plt.title('hyperopt first step')\n", + "plt.title(r'Restrict $s$ with polynomial')\n", "plt.legend()\n", "print('The minimum for cost function in the tested range is:', step_grid)" ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "Hence, we see that the strategy is indeed effective for finding the first minimum of the loss funciton for both the Z operator and the ZI operator." + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Compare in Pauli-Z strategy" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from qibo.quantum_info import random_hermitian\n", + "from qibo.hamiltonians import Hamiltonian" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Hamiltonian\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", + "nqubits = 4\n", + "h0 = random_hermitian(2**nqubits)\n", + "\n", + "# initialize class\n", + "dbi = DoubleBracketIteration(deepcopy(Hamiltonian(nqubits=nqubits, matrix=h0)),mode=DoubleBracketGeneratorType.single_commutator)\n", + "print(\"Initial off diagonal norm\", dbi.off_diagonal_norm)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "generate_local_Z = generate_Z_operators(nqubits)\n", + "Z_ops = list(generate_local_Z.values())\n", + "Z_names = list(generate_local_Z.keys())" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NSTEPS = 3\n", + "scheduling_list = [DoubleBracketScheduling.grid_search,\n", + " DoubleBracketScheduling.hyperopt,\n", + " DoubleBracketScheduling.polynomial_approximation,\n", + " DoubleBracketScheduling.simulated_annealing,]\n", + "scheduling_labels = ['grid search',\n", + " 'hyperopt',\n", + " 'polynomial',\n", + " 'simulated_annealing']\n", + "Z_optimal_scheduling = []\n", + "s_scheduling = []\n", + "off_norm_scheduling =[]\n", + "for i,scheduling in enumerate(scheduling_list):\n", + " # reinitialize\n", + " dbi = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=deepcopy(h0)), mode=DoubleBracketGeneratorType.single_commutator)\n", + " Z_optimal = []\n", + " # add in initial values for plotting\n", + " off_diagonal_norm_history = [dbi.off_diagonal_norm]\n", + " steps = [0]\n", + " print(f'----------Scheduling {scheduling_labels[i]}----------')\n", + " for _ in range(NSTEPS):\n", + " dbi, idx, step, flip_sign = select_best_dbr_generator(dbi, Z_ops, scheduling=scheduling, compare_canonical=False)\n", + " off_diagonal_norm_history.append(dbi.off_diagonal_norm)\n", + " steps.append(steps[-1]+step)\n", + " if flip_sign < 0:\n", + " Z_optimal.append('-' + Z_names[idx])\n", + " else:\n", + " Z_optimal.append(Z_names[idx])\n", + " print(f\"New optimized step at iteration {_+1}/{NSTEPS}: {step} with operator {Z_optimal[-1]}, loss {dbi.off_diagonal_norm}\")\n", + " Z_optimal_scheduling.append(Z_optimal)\n", + " s_scheduling.append(steps)\n", + " off_norm_scheduling.append(off_diagonal_norm_history)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.figure()\n", + "for i, scheduling in enumerate(scheduling_labels):\n", + " plt.plot(s_scheduling[i], off_norm_scheduling[i], '-o', label=scheduling)\n", + "plt.xlabel(\"Step durations\")\n", + "plt.ylabel(\"Norm off-diagonal restriction\")\n", + "plt.title(\"Compare Variational Pauli-Z using different scheduling strategies\")\n", + "plt.legend()" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## When polynomial approximation has no solution\n", + "\n", + "In some cases, the prescribed taylor expansion order `n` may not be sufficient to produce a meaningful step duration (real positive). In these cases, we rely on a backup scheduling method in `choose_step`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Hamiltonian\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", + "\n", + "# hamiltonian parameters\n", + "nqubits = 5\n", + "h = 3\n", + "\n", + "# define the hamiltonian\n", + "H_TFIM = hamiltonians.TFIM(nqubits=nqubits, h=h)\n", + "\n", + "# initialize class\n", + "dbi = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.canonical)\n", + "dbi.scheduling = DoubleBracketScheduling.polynomial_approximation\n", + "print(\"Initial off diagonal norm\", dbi.off_diagonal_norm)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "For demonstration purposes, we let `n=1` which is a linear fit to the loss function. This results in no valid solutions and function `polynomial_step` returns `None`." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "for n in range (5):\n", + " step = polynomial_step(dbi, n=n)\n", + " print(n, step)\n", + "print(dbi.choose_step(n=1))" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] } ], - "metadata": { - "kernelspec": { - "display_name": "DBF_qibo", - "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.11.7" - } - }, + "metadata": {}, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/examples/dbi/dbi_strategies_compare.ipynb b/examples/dbi/dbi_strategies_compare.ipynb new file mode 100644 index 0000000000..4a60ea3034 --- /dev/null +++ b/examples/dbi/dbi_strategies_compare.ipynb @@ -0,0 +1,464 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# DBI strategies comparison\n", + "\n", + "This notebook is a comparison of the so-far developed diagonalization strategies for DBI, including the canonical, Pauli-Z, and magnetic field strategies. On top of these, we also show case the use of invariant DBI generators such as 'BHMM'." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "from copy import deepcopy\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "from qibo import symbols\n", + "from qibo import hamiltonians, set_backend\n", + "from qibo.hamiltonians import Hamiltonian, SymbolicHamiltonian\n", + "from qibo.quantum_info import random_hermitian\n", + "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketScheduling, DoubleBracketIteration, DoubleBracketCostFunction\n", + "from qibo.models.dbi.utils import generate_Z_operators, generate_pauli_operator_dict, decompose_into_pauli_basis, ParameterizationTypes\n", + "from qibo.models.dbi.utils_dbr_strategies import select_best_dbr_generator, gradient_descent, polynomial_step" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def visualize_matrix(matrix, title=\"\"):\n", + " \"\"\"Visualize hamiltonian in a heatmap form.\"\"\"\n", + " fig, ax = plt.subplots(figsize=(5,5))\n", + " ax.set_title(title)\n", + " try:\n", + " im = ax.imshow(np.absolute(matrix), cmap=\"inferno\")\n", + " except TypeError:\n", + " im = ax.imshow(np.absolute(matrix.get()), cmap=\"inferno\")\n", + " fig.colorbar(im, ax=ax)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test on random Hamiltonian\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# backend\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", + "# initialize dbi object\n", + "nqubits = 5\n", + "h0 = random_hermitian(2**nqubits, seed=2)\n", + "dbi = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0))\n", + "cost = DoubleBracketCostFunction.off_diagonal_norm\n", + "print(\"Initial loss\", dbi.least_squares(d=dbi.diagonal_h_matrix))\n", + "visualize_matrix(dbi.h.matrix, title=f'Random hamiltonian with L={nqubits}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# iterations steps\n", + "NSTEPS = 5\n", + "# choose polynomial scheduling\n", + "scheduling = DoubleBracketScheduling.simulated_annealing" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Canonical" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# initialize DBI class for the canonical case\n", + "dbi_canonical = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0), mode=DoubleBracketGeneratorType.canonical, scheduling=scheduling, cost=cost)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Canonical\n", + "loss_history_canonical = [dbi_canonical.off_diagonal_norm]\n", + "steps_canonical_plot = [0]\n", + "for s in range(NSTEPS):\n", + " # same settings as iteration from list\n", + " d = dbi.diagonal_h_matrix\n", + " step = dbi_canonical.choose_step(d=d)\n", + " dbi_canonical(step=step)\n", + " print(f\"New optimized step at iteration {s+1}/{NSTEPS}: {step}, loss {dbi_canonical.off_diagonal_norm}\")\n", + " loss_history_canonical.append(dbi_canonical.off_diagonal_norm)\n", + " steps_canonical_plot.append(steps_canonical_plot[-1]+step)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Pauli-Z" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# initialize DBI class for the Pauli-Z strategy\n", + "dbi_pauli = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0), mode=DoubleBracketGeneratorType.single_commutator, scheduling=scheduling, cost=cost)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "generate_local_Z = generate_Z_operators(nqubits)\n", + "Z_ops = list(generate_local_Z.values())\n", + "Z_names = list(generate_local_Z.keys())\n", + "Z_optimal = []\n", + "# add in initial values for plotting\n", + "loss_history_pauli = [dbi_pauli.off_diagonal_norm]\n", + "steps_pauli_plot = [0]\n", + "scheduling = DoubleBracketScheduling.simulated_annealing\n", + "for _ in range(NSTEPS):\n", + " dbi_pauli, idx, step, flip_sign = select_best_dbr_generator(dbi_pauli, Z_ops, scheduling=scheduling, compare_canonical=False)\n", + " d = Z_ops[idx]\n", + " loss_history_pauli.append(dbi_pauli.off_diagonal_norm)\n", + " steps_pauli_plot.append(steps_pauli_plot[-1]+step)\n", + " if flip_sign < 0:\n", + " Z_optimal.append('-' + Z_names[idx])\n", + " else:\n", + " Z_optimal.append(Z_names[idx])\n", + " print(f\"New optimized step at iteration {_+1}/{NSTEPS}: {step} with operator {Z_optimal[-1]}, loss {dbi_pauli.off_diagonal_norm}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Magnetic field" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# initialize DBI class for the canonical case\n", + "dbi_gradient = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0), mode=DoubleBracketGeneratorType.single_commutator, scheduling=scheduling, cost=cost)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pauli_operator_dict = generate_pauli_operator_dict(nqubits=nqubits, parameterization_order=2)\n", + "d_coef = decompose_into_pauli_basis(dbi.h.matrix, list(pauli_operator_dict.values()))\n", + "d = sum([d_coef[i] * list(pauli_operator_dict.values())[i] for i in range(nqubits)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def s_hist_to_plot(s_hist):\n", + " # convert list of step durations taken to plotable\n", + " s_plot = [0] * len(s_hist)\n", + " for i in range(len(s_hist)):\n", + " if i != 0:\n", + " s_plot[i] = s_plot[i-1] + s_hist[i]\n", + " return s_plot" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "loss_history_gradient, d_params_hist, s_hist = gradient_descent(dbi_gradient, NSTEPS, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)\n", + "steps_gradient_plot = s_hist_to_plot(s_hist)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(str(nqubits) + ' random Hamiltonian diagonalization')\n", + "plt.plot(loss_history_canonical, label='canonical')\n", + "plt.plot(loss_history_pauli, label='Pauli-Z')\n", + "plt.plot(loss_history_gradient, label='gradient')\n", + "plt.legend()\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(str(nqubits) + ' random Hamiltonian diagonalization')\n", + "plt.plot(steps_canonical_plot, loss_history_canonical, marker='o', label='canonical')\n", + "plt.plot(steps_pauli_plot, loss_history_pauli, marker='o', label='Pauli-Z')\n", + "plt.plot(steps_gradient_plot,loss_history_gradient, marker='o', label='gradient')\n", + "plt.legend()\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test on TFIM\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# backend\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", + "# initialize dbi object\n", + "# hamiltonian parameters\n", + "# define the hamiltonian\n", + "nqubits = 5\n", + "h = 1\n", + "H_TFIM = SymbolicHamiltonian( - h*symbols.Z(nqubits-1), nqubits=nqubits)\n", + "# add linear interaction terms\n", + "for i in range(nqubits-1):\n", + " H_TFIM -= SymbolicHamiltonian(symbols.X(i)*symbols.X(i+1) + h*symbols.Z(i), nqubits=nqubits)\n", + "H_TFIM = H_TFIM.dense\n", + "\n", + "# initialize class\n", + "dbi = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.canonical)\n", + "print(\"Initial off diagonal norm\", dbi.off_diagonal_norm)\n", + "visualize_matrix(dbi.h.matrix, title=f'Random hamiltonian with L={nqubits}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# iterations steps\n", + "NSTEPS = 5\n", + "# choose polynomial scheduling\n", + "scheduling = DoubleBracketScheduling.simulated_annealing" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Canonical" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# initialize DBI class for the canonical case\n", + "dbi_canonical = DoubleBracketIteration(deepcopy(H_TFIM), mode=DoubleBracketGeneratorType.canonical, scheduling=scheduling)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# Canonical\n", + "off_diagonal_norm_history_canonical = [dbi_canonical.off_diagonal_norm]\n", + "steps_canonical_plot = [0]\n", + "for s in range(NSTEPS):\n", + " # same settings as iteration from list\n", + " step = dbi_canonical.choose_step(d=dbi.diagonal_h_matrix)\n", + " dbi_canonical(step=step)\n", + " print(f\"New optimized step at iteration {s+1}/{NSTEPS}: {step}, loss {dbi_canonical.off_diagonal_norm}\")\n", + " off_diagonal_norm_history_canonical.append(dbi_canonical.off_diagonal_norm)\n", + " steps_canonical_plot.append(steps_canonical_plot[-1]+step)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Pauli-Z" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# initialize DBI class for the Pauli-Z strategy\n", + "dbi_pauli = DoubleBracketIteration(deepcopy(H_TFIM), mode=DoubleBracketGeneratorType.single_commutator, scheduling=scheduling)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "generate_local_Z = generate_Z_operators(nqubits)\n", + "Z_ops = list(generate_local_Z.values())\n", + "Z_names = list(generate_local_Z.keys())\n", + "Z_optimal = []\n", + "# add in initial values for plotting\n", + "off_diagonal_norm_history_pauli = [dbi_pauli.off_diagonal_norm]\n", + "steps_pauli_plot = [0]\n", + "scheduling = DoubleBracketScheduling.simulated_annealing\n", + "for _ in range(NSTEPS):\n", + " dbi_pauli, idx, step, flip_sign = select_best_dbr_generator(dbi_pauli, Z_ops, scheduling=scheduling, compare_canonical=False)\n", + " off_diagonal_norm_history_pauli.append(dbi_pauli.off_diagonal_norm)\n", + " steps_pauli_plot.append(steps_pauli_plot[-1]+step)\n", + " if flip_sign < 0:\n", + " Z_optimal.append('-' + Z_names[idx])\n", + " else:\n", + " Z_optimal.append(Z_names[idx])\n", + " print(f\"New optimized step at iteration {_+1}/{NSTEPS}: {step} with operator {Z_optimal[-1]}, loss {dbi_pauli.off_diagonal_norm}\")" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Magnetic field" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# initialize DBI class for the canonical case\n", + "dbi_gradient = DoubleBracketIteration(deepcopy(H_TFIM), mode=DoubleBracketGeneratorType.single_commutator, scheduling=scheduling)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "pauli_operator_dict = generate_pauli_operator_dict(nqubits=nqubits, parameterization_order=2)\n", + "d_coef = decompose_into_pauli_basis(dbi.h.matrix, list(pauli_operator_dict.values()))\n", + "d = sum([d_coef[i] * list(pauli_operator_dict.values())[i] for i in range(nqubits)])" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "initial_s = polynomial_step(dbi_object=dbi, d=d, n=4)\n", + "print(initial_s)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "off_diagonal_norm_history_gradient, d_params_hist, s_hist = gradient_descent(dbi_gradient, NSTEPS, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)\n", + "steps_gradient_plot = s_hist_to_plot(s_hist)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(str(nqubits) + ' random Hamiltonian diagonalization')\n", + "plt.plot(off_diagonal_norm_history_canonical, label='canonical', marker='o')\n", + "plt.plot(off_diagonal_norm_history_pauli, label='Pauli-Z', marker='o')\n", + "plt.plot(off_diagonal_norm_history_gradient, label='gradient', marker='o')\n", + "plt.legend()\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(str(nqubits) + ' random Hamiltonian diagonalization')\n", + "plt.plot(steps_canonical_plot, off_diagonal_norm_history_canonical, marker='o', label='canonical')\n", + "plt.plot(steps_pauli_plot, off_diagonal_norm_history_pauli, marker='o', label='Pauli-Z')\n", + "plt.plot(steps_gradient_plot,off_diagonal_norm_history_gradient, marker='o', label='gradient')\n", + "plt.legend()\n", + "plt.xlabel('Duration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/dbi/dbi_strategy_Ising_model.ipynb b/examples/dbi/dbi_strategy_Ising_model.ipynb new file mode 100644 index 0000000000..6c0d53209a --- /dev/null +++ b/examples/dbi/dbi_strategy_Ising_model.ipynb @@ -0,0 +1,328 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "# Double-Bracket Iteration Strategy: magnetic field (Ising model)\n", + "This notebook shows the diagonalization process of DBI using the magnetic field strategy, which varies the diagonal operator $D$ by gradient descent. To find the gradient with respect to $D$, parameterization of $D$ is required. For the purpose of this notebook, we represent it by the Ising model, i.e.\n", + "\n", + "$$ D = \\sum \\alpha_i Z_i +\\sum \\beta_{ij}Z_iZ_j$$\n", + "\n", + "\n", + "The gradients are calculated under the premise that the diagonalization gain curve can be fitted by a polynomial, and that the iteration step duration is taken at the first dip of the curve." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "from copy import deepcopy\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "from qibo import hamiltonians, set_backend, symbols\n", + "from qibo.hamiltonians import Hamiltonian, SymbolicHamiltonian\n", + "from qibo.quantum_info import random_hermitian\n", + "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketScheduling, DoubleBracketIteration\n", + "from qibo.models.dbi.utils import generate_pauli_operator_dict, decompose_into_pauli_basis, ParameterizationTypes\n", + "from qibo.models.dbi.utils_dbr_strategies import gradient_numerical, gradient_descent" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def visualize_matrix(matrix, title=\"\"):\n", + " \"\"\"Visualize hamiltonian in a heatmap form.\"\"\"\n", + " fig, ax = plt.subplots(figsize=(5,5))\n", + " ax.set_title(title)\n", + " try:\n", + " im = ax.imshow(np.absolute(matrix), cmap=\"inferno\")\n", + " except TypeError:\n", + " im = ax.imshow(np.absolute(matrix.get()), cmap=\"inferno\")\n", + " fig.colorbar(im, ax=ax)\n" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test on random Hamiltonian" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# backend\n", + "set_backend(\"numpy\")\n", + "# initialize dbi object\n", + "nqubits = 5\n", + "h0 = random_hermitian(2**nqubits, seed=2)\n", + "scheduling = DoubleBracketScheduling.hyperopt\n", + "mode = DoubleBracketGeneratorType.single_commutator\n", + "n_taylor = 5\n", + "dbi = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0), scheduling=scheduling, mode=mode)\n", + "print(\"Initial off diagonal norm\", dbi.off_diagonal_norm)\n", + "visualize_matrix(dbi.h.matrix, title=f'Random hamiltonian with L={nqubits}')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Order 1: $D=\\sum \\alpha_iZ_i$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# generate pauli_operator_dict\n", + "pauli_operator_dict = generate_pauli_operator_dict(nqubits=nqubits, parameterization_order=1)\n", + "d_coef = decompose_into_pauli_basis(dbi.h.matrix, list(pauli_operator_dict.values()))\n", + "d = sum([d_coef[i] * list(pauli_operator_dict.values())[i] for i in range(nqubits)])\n", + "grad = gradient_numerical(dbi, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)\n", + "print('The initial D coefficients:', d_coef)\n", + "print('Gradient:', grad)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iters = 15\n", + "off_diagonal_norm_1, d_params_hist, s_step = gradient_descent(dbi, iters, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(str(nqubits) + ' spins random hamiltonian')\n", + "plt.plot(off_diagonal_norm_1)\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Order 2: $D=\\sum \\alpha_iZ_i + \\beta_{ij}Z_iZ_j$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dbi = DoubleBracketIteration(Hamiltonian(nqubits=nqubits, matrix=h0), scheduling=scheduling, mode=mode)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# generate pauli_operator_dict\n", + "pauli_operator_dict = generate_pauli_operator_dict(nqubits=nqubits, parameterization_order=2)\n", + "d_coef = decompose_into_pauli_basis(dbi.h.matrix, list(pauli_operator_dict.values()))\n", + "d = sum([d_coef[i] * list(pauli_operator_dict.values())[i] for i in range(nqubits)])\n", + "grad = gradient_numerical(dbi, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)\n", + "print('The initial D coefficients:', d_coef)\n", + "print('Gradient:', grad)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "iters = 15\n", + "off_diagonal_norm_2, d_params_hist, s_step = gradient_descent(dbi, iters, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(str(nqubits) + ' spins random hamiltonian')\n", + "plt.plot(off_diagonal_norm_1, label='order 1')\n", + "plt.plot(off_diagonal_norm_2, label='order 2')\n", + "plt.legend()\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "## Test on TFIM\n", + "Here we choose to customize our TFIM in the X axis using `SymbolicHamiltonian`. It is also possible to use Hadamard gate to rotate the TFIM inbuilt in `qibo`.\n", + "\n", + "$$ H = -(\\sum X_i X_{i+1} + \\sum hZ_i)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# generate the Hamiltonian\n", + "nqubits = 5\n", + "h = 1\n", + "H_TFIM = SymbolicHamiltonian( - h*symbols.Z(nqubits-1), nqubits=nqubits)\n", + "# add linear interaction terms\n", + "for i in range(nqubits-1):\n", + " H_TFIM -= SymbolicHamiltonian(symbols.X(i)*symbols.X(i+1) + h*symbols.Z(i), nqubits=nqubits)\n", + "H_TFIM = H_TFIM.dense\n", + "visualize_matrix(H_TFIM.matrix, title=f'TFIM with L={nqubits} h={h}')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# backend\n", + "set_backend(\"numpy\")\n", + "# initialize dbi object\n", + "dbi_TFIM = DoubleBracketIteration(deepcopy(H_TFIM), scheduling=scheduling, mode=mode)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Order 1: $D=\\sum \\alpha_iZ_i$" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dbi_TFIM_1 = DoubleBracketIteration(deepcopy(H_TFIM), scheduling=scheduling, mode=mode)\n", + "# generate pauli_operator_dict\n", + "pauli_operator_dict = generate_pauli_operator_dict(nqubits=nqubits, parameterization_order=1)\n", + "d_coef = decompose_into_pauli_basis(dbi_TFIM_1.h.matrix, list(pauli_operator_dict.values()))\n", + "d = sum([d_coef[i] * list(pauli_operator_dict.values())[i] for i in range(nqubits)])\n", + "grad = gradient_numerical(dbi_TFIM_1, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)\n", + "print('The initial D coefficients:', d_coef)\n", + "print('Gradient:', grad)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "NSTEPS = 3\n", + "off_diagonal_norm_1, d_params_hist, s_step = gradient_descent(dbi_TFIM_1, NSTEPS, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(f'n={nqubits} h={h} TFIM, order=1')\n", + "plt.plot(off_diagonal_norm_1)\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# the final matrix\n", + "visualize_matrix(dbi_TFIM.h.matrix)" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "### Order 2" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "dbi_TFIM_2 = DoubleBracketIteration(deepcopy(H_TFIM), scheduling=scheduling, mode=mode)\n", + "# generate pauli_operator_dict\n", + "pauli_operator_dict = generate_pauli_operator_dict(nqubits=nqubits, parameterization_order=2)\n", + "d_coef = decompose_into_pauli_basis(dbi_TFIM_2.h.matrix, list(pauli_operator_dict.values()))\n", + "d = sum([d_coef[i] * list(pauli_operator_dict.values())[i] for i in range(nqubits)])\n", + "grad = gradient_numerical(dbi_TFIM_2, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)\n", + "print('The initial D coefficients:', d_coef)\n", + "print('Gradient:', grad)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "off_diagonal_norm_2, d_params_hist, s_step = gradient_descent(dbi_TFIM_2, NSTEPS, d_coef, parameterization=ParameterizationTypes.pauli, pauli_operator_dict=pauli_operator_dict)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "plt.title(f'n={nqubits} h={h} TFIM')\n", + "plt.plot(off_diagonal_norm_1, label='order 1')\n", + "plt.plot(off_diagonal_norm_2, label='order 2')\n", + "plt.legend()\n", + "plt.xlabel('Iteration')\n", + "plt.ylabel(r'$|| \\sigma(e^{sW}He^{-sW}) || $')" + ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [ + "In conclusion, we see that the parameterization order or locality of the Pauli based parameterization for gradient descent does not affect significantly the effectiveness of double bracket diagonalization." + ] + } + ], + "metadata": {}, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/dbi/dbi_strategy_Pauli-Z.ipynb b/examples/dbi/dbi_strategy_Pauli-Z.ipynb index 0f76a36245..aabee6c172 100644 --- a/examples/dbi/dbi_strategy_Pauli-Z.ipynb +++ b/examples/dbi/dbi_strategy_Pauli-Z.ipynb @@ -22,13 +22,12 @@ "metadata": {}, "outputs": [], "source": [ - "!python -m pip install seaborn # plotting library\n", - "!python -m pip install hyperopt # required to optimize the DBF step" + "# !python -m pip install seaborn # plotting library" ] }, { "cell_type": "code", - "execution_count": 2, + "execution_count": null, "metadata": {}, "outputs": [], "source": [ @@ -39,8 +38,9 @@ "import seaborn as sns\n", "\n", "from qibo import hamiltonians, set_backend\n", - "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketIteration\n", - "from qibo.models.dbi.utils import *" + "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketIteration, DoubleBracketScheduling\n", + "from qibo.models.dbi.utils import generate_Z_operators\n", + "from qibo.models.dbi.utils_dbr_strategies import select_best_dbr_generator" ] }, { @@ -122,31 +122,16 @@ }, { "cell_type": "code", - "execution_count": 9, - "metadata": {}, - "outputs": [ - { - "name": "stderr", - "output_type": "stream", - "text": [ - "[Qibo 0.2.4|INFO|2024-01-24 19:59:31]: Using qibojit (numba) backend on /CPU:0\n" - ] - }, - { - "name": "stdout", - "output_type": "stream", - "text": [ - "Initial off diagonal norm 8.48528137423857\n" - ] - } - ], + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "# set the qibo backend (we suggest qibojit if N >= 20)\n", "# alternatives: tensorflow (not optimized), numpy (when CPU not supported by jit)\n", - "set_backend(\"qibojit\", \"numba\")\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", "\n", "# hamiltonian parameters\n", - "nqubits = 2\n", + "nqubits = 5\n", "h = 3\n", "\n", "# define the hamiltonian\n", @@ -160,20 +145,9 @@ }, { "cell_type": "code", - "execution_count": 8, - "metadata": {}, - "outputs": [ - { - "name": "stdout", - "output_type": "stream", - "text": [ - "[[-2.-0.j -0.-0.j -0.-0.j -0.-0.j]\n", - " [-0.-0.j 2.-0.j -0.-0.j -0.-0.j]\n", - " [-0.-0.j -0.-0.j 2.-0.j -0.-0.j]\n", - " [-0.-0.j -0.-0.j -0.-0.j -2.-0.j]]\n" - ] - } - ], + "execution_count": null, + "metadata": {}, + "outputs": [], "source": [ "print(H_TFIM.matrix)" ] @@ -212,15 +186,16 @@ "metadata": {}, "outputs": [], "source": [ - "NSTEPS = 15\n", - "max_evals = 100\n", + "NSTEPS = 2\n", + "max_evals = 10\n", "step_max = 1\n", "Z_optimal = []\n", "# add in initial values for plotting\n", "off_diagonal_norm_history = [dbi.off_diagonal_norm]\n", "steps = [0]\n", + "scheduling = DoubleBracketScheduling.hyperopt\n", "for _ in range(NSTEPS):\n", - " dbi, idx, step, flip_sign = select_best_dbr_generator(dbi, Z_ops, compare_canonical=False, max_evals=max_evals, step_max=step_max)\n", + " dbi, idx, step, flip_sign = select_best_dbr_generator(dbi, Z_ops, scheduling=scheduling, compare_canonical=False, max_evals=max_evals, step_max=step_max)\n", " off_diagonal_norm_history.append(dbi.off_diagonal_norm)\n", " steps.append(steps[-1]+step)\n", " if flip_sign < 0:\n", @@ -269,7 +244,7 @@ "source": [ "# set the qibo backend (we suggest qibojit if N >= 20)\n", "# alternatives: tensorflow (not optimized), numpy (when CPU not supported by jit)\n", - "set_backend(\"qibojit\", \"numba\")\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", "\n", "\n", "# initialize class|\n", @@ -289,13 +264,7 @@ "steps_canonical_plot = [0]\n", "for s in range(NSTEPS):\n", " # same settings as iteration from list\n", - " step = dbi_canonical.hyperopt_step(\n", - " step_min = 1e-5,\n", - " step_max = 1,\n", - " space = hp.uniform,\n", - " optimizer = tpe,\n", - " max_evals = max_evals,\n", - " )\n", + " step = dbi_canonical.choose_step(scheduling=DoubleBracketScheduling.hyperopt)\n", " dbi_canonical(step=step)\n", " print(f\"New optimized step at iteration {s+1}/{NSTEPS}: {step}, loss {dbi_canonical.off_diagonal_norm}\")\n", " off_diagonal_norm_history_canonical.append(dbi_canonical.off_diagonal_norm)\n", @@ -355,7 +324,7 @@ "metadata": {}, "outputs": [], "source": [ - "dbi_mixed = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.single_commutator)\n", + "dbi_mixed = DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.single_commutator, scheduling=DoubleBracketScheduling.hyperopt)\n", "print(\"Initial off diagonal norm\", dbi_mixed.off_diagonal_norm)" ] }, @@ -368,12 +337,7 @@ "dbi_eval = deepcopy(dbi_mixed)\n", "dbi_eval.mode = DoubleBracketGeneratorType.canonical\n", "if step is None:\n", - " step = dbi_eval.hyperopt_step(\n", - " step_max=step_max,\n", - " space=hp.uniform,\n", - " optimizer=tpe,\n", - " max_evals=max_evals,\n", - " )\n", + " step = dbi_eval.choose_step()\n", "dbi_eval(step=step)\n", "print('canonical norm', dbi_eval.off_diagonal_norm, 'step', step)" ] @@ -389,7 +353,7 @@ "off_diagonal_norm_history_mixed = [dbi_mixed.off_diagonal_norm]\n", "steps = [0]\n", "for _ in range(NSTEPS):\n", - " dbi_mixed, idx, step, flip_sign = select_best_dbr_generator(dbi_mixed, Z_ops, compare_canonical=True, max_evals=max_evals)\n", + " dbi_mixed, idx, step, flip_sign = select_best_dbr_generator(dbi_mixed, Z_ops, scheduling=scheduling, compare_canonical=True, max_evals=max_evals, step_max=step_max)\n", " off_diagonal_norm_history_mixed.append(dbi_mixed.off_diagonal_norm)\n", " steps.append(steps[-1]+step)\n", " if idx == len(Z_ops):\n", @@ -444,7 +408,7 @@ "metadata": {}, "outputs": [], "source": [ - "dbi_mixed_can= DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.canonical)\n", + "dbi_mixed_can= DoubleBracketIteration(deepcopy(H_TFIM),mode=DoubleBracketGeneratorType.canonical, scheduling=DoubleBracketScheduling.hyperopt)\n", "print(\"Initial off diagonal norm\", dbi_mixed_can.off_diagonal_norm)" ] }, @@ -463,7 +427,7 @@ " dbi_mixed_can(step=step)\n", " off_diagonal_norm_history_mixed_can.append(dbi_mixed_can.off_diagonal_norm)\n", " steps_mixed_can.append(step)\n", - " \n", + "\n", "print(\"After 2 steps, off diagonal norm:\", dbi_mixed_can.off_diagonal_norm)\n", "print(\"By comparison, the Pauli-Z:\", off_diagonal_norm_history[2])" ] @@ -479,7 +443,7 @@ "remaining_NSTEPS = NSTEPS - cannonical_NSTEPS\n", "dbi_mixed_can.mode = DoubleBracketGeneratorType.single_commutator\n", "for _ in range(remaining_NSTEPS):\n", - " dbi_mixed_can, idx, step, flip_sign = select_best_dbr_generator(dbi_mixed_can, Z_ops, compare_canonical=False)\n", + " dbi_mixed_can, idx, step, flip_sign = select_best_dbr_generator(dbi_mixed_can, Z_ops, scheduling=scheduling, compare_canonical=False, max_evals=max_evals, step_max=step_max)\n", " off_diagonal_norm_history_mixed_can.append(dbi_mixed_can.off_diagonal_norm)\n", " steps_mixed_can.append(step)\n", " if idx == len(Z_ops):\n", @@ -516,27 +480,14 @@ "source": [ "This example also shows that the canonical generator is more likely to drive the model into a local minimum than variationally assigned diagonal operator, and that it is hard to get it unstuck even with the Pauli-Z operators." ] + }, + { + "cell_type": "markdown", + "metadata": {}, + "source": [] } ], - "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.11.7" - } - }, + "metadata": {}, "nbformat": 4, - "nbformat_minor": 2 + "nbformat_minor": 4 } diff --git a/examples/dbi/dbi_tutorial_basic_intro.ipynb b/examples/dbi/dbi_tutorial_basic_intro.ipynb index ecb28bb4d7..bcd0d65c67 100644 --- a/examples/dbi/dbi_tutorial_basic_intro.ipynb +++ b/examples/dbi/dbi_tutorial_basic_intro.ipynb @@ -2,7 +2,7 @@ "cells": [ { "cell_type": "markdown", - "id": "cb748c1a-2ecd-44a2-91d8-c1255a00615b", + "id": "2a33581d", "metadata": {}, "source": [ "## Double-Bracket Iteration diagonalization algorithm\n", @@ -17,18 +17,18 @@ { "cell_type": "code", "execution_count": null, - "id": "e1f362b8-eb73-456e-ae48-94c5f2a12649", + "id": "62d9723f", "metadata": {}, "outputs": [], "source": [ - "!python -m pip install seaborn # plotting library\n", - "!python -m pip install hyperopt # required to optimize the DBF step" + "# uncomment this line if seaborn is not installed\n", + "# !python -m pip install seaborn" ] }, { "cell_type": "code", "execution_count": null, - "id": "f270b1ea-ee6a-4eac-a0ff-3d7dae296cf0", + "id": "b80b4738", "metadata": {}, "outputs": [], "source": [ @@ -41,12 +41,12 @@ "from hyperopt import hp, tpe\n", "\n", "from qibo import hamiltonians, set_backend\n", - "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketIteration" + "from qibo.models.dbi.double_bracket import DoubleBracketGeneratorType, DoubleBracketIteration, DoubleBracketScheduling" ] }, { "cell_type": "markdown", - "id": "ba6e5402-ea34-4979-bb79-fd395567f77d", + "id": "a5e25f51", "metadata": {}, "source": [ "Here we define a simple plotting function useful to keep track of the diagonalization process." @@ -55,7 +55,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4aec7b46-19b9-4004-93c0-a90255e58fd9", + "id": "933d9a00", "metadata": {}, "outputs": [], "source": [ @@ -96,7 +96,7 @@ }, { "cell_type": "markdown", - "id": "9f4cd7cc-9952-4da4-baef-e916300a9365", + "id": "4efd4a97", "metadata": {}, "source": [ "We need to define a target hamiltonian which we aim to diagonalize. As an example, we consider the Transverse Field Ising Model (TFIM):\n", @@ -107,12 +107,12 @@ { "cell_type": "code", "execution_count": null, - "id": "2c4ed408-68ed-4054-825c-2a7df0979a4f", + "id": "7125940f", "metadata": {}, "outputs": [], "source": [ "# set the qibo backend (we suggest qibojit if N >= 20)\n", - "set_backend(\"qibojit\", \"numba\")\n", + "set_backend(\"qibojit\", platform=\"numba\")\n", "\n", "# hamiltonian parameters\n", "nqubits = 5\n", @@ -127,7 +127,7 @@ }, { "cell_type": "markdown", - "id": "4794e779-bf2d-4ab5-97ce-f876d9522a35", + "id": "c2ca8392", "metadata": {}, "source": [ "#### The generator of the evolution\n", @@ -148,7 +148,7 @@ { "cell_type": "code", "execution_count": null, - "id": "26a487e9-366b-4203-b660-e3d4af2bcb68", + "id": "1adafc19", "metadata": {}, "outputs": [], "source": [ @@ -160,7 +160,7 @@ { "cell_type": "code", "execution_count": null, - "id": "da8dce89-27f6-403d-982a-58d531fade48", + "id": "8a4d0e9d", "metadata": {}, "outputs": [], "source": [ @@ -170,7 +170,7 @@ }, { "cell_type": "markdown", - "id": "fc4f9f75-0548-4533-a13c-3aba3191e608", + "id": "a5527622", "metadata": {}, "source": [ "#### The `DoubleBracketIteration` class\n", @@ -181,7 +181,7 @@ { "cell_type": "code", "execution_count": null, - "id": "055870ec-55f2-4b99-a622-e3aa4c7dd0e9", + "id": "9521c464", "metadata": {}, "outputs": [], "source": [ @@ -190,7 +190,7 @@ }, { "cell_type": "markdown", - "id": "b38cf803-60b4-467a-be8e-cbad5d81f14a", + "id": "a262c69f", "metadata": {}, "source": [ "#### `DoubleBracketIteration` features" @@ -199,7 +199,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9e278c3d-9f34-4a40-b453-4e030c751ef5", + "id": "290e5828", "metadata": {}, "outputs": [], "source": [ @@ -210,7 +210,7 @@ { "cell_type": "code", "execution_count": null, - "id": "5b8e142b-a0a2-41bd-a16a-265a420b7360", + "id": "3e2b9950", "metadata": {}, "outputs": [], "source": [ @@ -222,7 +222,7 @@ { "cell_type": "code", "execution_count": null, - "id": "4f9d1d41-3df7-49cf-96ca-fa1019c00c33", + "id": "638ba4b5", "metadata": {}, "outputs": [], "source": [ @@ -233,7 +233,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7b864712-219c-44b6-8337-19ef0100e318", + "id": "08f0c466", "metadata": {}, "outputs": [], "source": [ @@ -243,7 +243,7 @@ }, { "cell_type": "markdown", - "id": "5e576bc4-4e79-4c71-9ea0-b3012e9f2ba1", + "id": "bb5f10da", "metadata": {}, "source": [ "which shows $\\hat{H}$ is now identical to $\\hat{H}_0$ since no evolution step has been performed yet." @@ -252,7 +252,7 @@ { "cell_type": "code", "execution_count": null, - "id": "da3d3aaa-17e1-492e-bcd3-b510f44a5391", + "id": "90e6fdff", "metadata": {}, "outputs": [], "source": [ @@ -262,7 +262,7 @@ }, { "cell_type": "markdown", - "id": "ca0ce252", + "id": "a0101ae0", "metadata": {}, "source": [ "The Hilbert-Schmidt norm of a Hamiltonian is defined as:\n", @@ -273,7 +273,7 @@ { "cell_type": "code", "execution_count": null, - "id": "24d0dfa1-7039-4d7d-8aa3-5a937b9ab0b8", + "id": "0d90c8b5", "metadata": {}, "outputs": [], "source": [ @@ -284,7 +284,7 @@ }, { "cell_type": "markdown", - "id": "d75e35ab-66f4-49f9-af19-679c20065a11", + "id": "a1d1eb77", "metadata": {}, "source": [ "Finally, the energy fluctuation of the system at step $k$ over a given state $\\mu$\n", @@ -297,7 +297,7 @@ { "cell_type": "code", "execution_count": null, - "id": "95f8d86f-07d4-498c-acb1-f6f6a4614c24", + "id": "13710cc2", "metadata": {}, "outputs": [], "source": [ @@ -312,7 +312,7 @@ }, { "cell_type": "markdown", - "id": "3d5b37f3-2477-49a0-9f80-7da5ddda1fff", + "id": "4d34e1e3", "metadata": {}, "source": [ "#### Call the `DoubleBracketIteration` to perform a DBF iteration\n", @@ -323,7 +323,7 @@ { "cell_type": "code", "execution_count": null, - "id": "9a886261-8aa6-4df0-a31b-9c39847db124", + "id": "a7749a96", "metadata": {}, "outputs": [], "source": [ @@ -340,7 +340,7 @@ }, { "cell_type": "markdown", - "id": "b78dd05d-ffe3-435a-b5ec-2a42f28066b2", + "id": "dab441bb", "metadata": {}, "source": [ "We can check now if something happened by plotting the drift:" @@ -349,7 +349,7 @@ { "cell_type": "code", "execution_count": null, - "id": "cc74812d-7c2c-44e4-afc2-e235968801b4", + "id": "fc01baa4", "metadata": {}, "outputs": [], "source": [ @@ -358,7 +358,7 @@ }, { "cell_type": "markdown", - "id": "3465a422-eebf-4e80-ae96-bba894132330", + "id": "9223433b", "metadata": {}, "source": [ "The set step can be good, but maybe not the best one. In order to do this choice in a wiser way, we can call the DBF hyperoptimization routine to search for a better initial step. The `dbf.hyperopt_step` method is built on top of the [`hyperopt`](https://hyperopt.github.io/hyperopt/) package. Any algorithm or sampling space provided by the official package can be used. We are going to use the default options (we sample new steps from a uniform space following a _Tree of Parzen estimators algorithm_)." @@ -367,7 +367,7 @@ { "cell_type": "code", "execution_count": null, - "id": "aad79966-7a11-4a45-aba5-4a4bb8315c50", + "id": "0d7b86d3", "metadata": {}, "outputs": [], "source": [ @@ -375,20 +375,20 @@ "dbf.h = dbf.h0\n", "\n", "# optimization of the step, we allow to search in [1e-5, 1]\n", - "step = dbf.hyperopt_step(\n", + "step = dbf.choose_step(\n", + " scheduling=DoubleBracketScheduling.hyperopt,\n", " step_min = 1e-5,\n", " step_max = 1,\n", " space = hp.uniform,\n", " optimizer = tpe,\n", " max_evals = 1000,\n", - " verbose = True\n", ")" ] }, { "cell_type": "code", "execution_count": null, - "id": "49483a47-d29d-440e-a4bc-143bfe6bb3cf", + "id": "1b9b1431", "metadata": {}, "outputs": [], "source": [ @@ -398,7 +398,7 @@ { "cell_type": "code", "execution_count": null, - "id": "6bdaf7f9-7e49-4a16-8b29-ae1f9746cd9b", + "id": "52fa3599", "metadata": {}, "outputs": [], "source": [ @@ -407,7 +407,7 @@ }, { "cell_type": "markdown", - "id": "b5f1d00e-e763-40d9-822f-e0e8d4c57d9a", + "id": "084c3bcb", "metadata": {}, "source": [ "#### Let's evolve the model for `NSTEPS`\n", @@ -420,7 +420,7 @@ { "cell_type": "code", "execution_count": null, - "id": "59a6a485-a714-4e14-b27a-1df2930068ee", + "id": "d1f197b1", "metadata": {}, "outputs": [], "source": [ @@ -443,7 +443,7 @@ { "cell_type": "code", "execution_count": null, - "id": "7e0b2f18-ca53-4f34-9fcf-0052dcc31dc5", + "id": "c115c222", "metadata": {}, "outputs": [], "source": [ @@ -452,7 +452,7 @@ }, { "cell_type": "markdown", - "id": "eb797d6c-0eba-4da4-b492-8b5d70f9123f", + "id": "233ba431", "metadata": {}, "source": [ "#### Method 2: optimizing the step" @@ -461,30 +461,29 @@ { "cell_type": "code", "execution_count": null, - "id": "a6fd1e33-3620-4f3b-b705-a120f6da0027", + "id": "4e0fc1c2", "metadata": {}, "outputs": [], "source": [ "# restart\n", - "dbf_2 = DoubleBracketIteration(hamiltonian=deepcopy(h), mode=iterationtype)\n", + "dbf_2 = DoubleBracketIteration(hamiltonian=deepcopy(h), mode=iterationtype, scheduling=DoubleBracketScheduling.hyperopt)\n", "off_diagonal_norm_history = [dbf_2.off_diagonal_norm]\n", "\n", "# set the number of evolution steps\n", "NSTEPS = 20\n", "\n", "# optimize first step\n", - "step = dbf_2.hyperopt_step(\n", + "step = dbf_2.choose_step(\n", " step_min = 1e-5,\n", " step_max = 1,\n", " space = hp.uniform,\n", " optimizer = tpe,\n", " max_evals = 500,\n", - " verbose = True\n", ")\n", "\n", "for s in range(NSTEPS):\n", " if s != 0:\n", - " step = dbf_2.hyperopt_step(\n", + " step = dbf_2.choose_step(\n", " step_min = 1e-5,\n", " step_max = 1,\n", " space = hp.uniform,\n", @@ -502,7 +501,7 @@ { "cell_type": "code", "execution_count": null, - "id": "0f0212bf-b642-4fea-9203-037876e0b266", + "id": "40e31e97", "metadata": {}, "outputs": [], "source": [ @@ -511,7 +510,7 @@ }, { "cell_type": "markdown", - "id": "32341937-4178-41d2-a10e-5e4d2634098e", + "id": "0de78acd", "metadata": {}, "source": [ "The hyperoptimization can lead to a faster convergence of the algorithm." @@ -520,7 +519,7 @@ { "cell_type": "code", "execution_count": null, - "id": "82b89092-07e5-4788-9ae0-8907df2428eb", + "id": "baab0ab5", "metadata": {}, "outputs": [], "source": [ @@ -530,33 +529,23 @@ { "cell_type": "code", "execution_count": null, - "id": "ac8ed320-04a8-42af-a980-48ab4f1fff7c", + "id": "2bc9ac69", "metadata": {}, "outputs": [], "source": [ "visualize_matrix(dbf_2.h.matrix)" ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0bed191d", + "metadata": {}, + "outputs": [], + "source": [] } ], - "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.11.7" - } - }, + "metadata": {}, "nbformat": 4, "nbformat_minor": 5 } diff --git a/poetry.lock b/poetry.lock index 37c33d53c9..dd45df2f26 100644 --- a/poetry.lock +++ b/poetry.lock @@ -55,6 +55,17 @@ doc = ["Sphinx (>=7)", "packaging", "sphinx-autodoc-typehints (>=1.2.0)", "sphin test = ["anyio[trio]", "coverage[toml] (>=7)", "exceptiongroup (>=1.2.0)", "hypothesis (>=4.0)", "psutil (>=5.9)", "pytest (>=7.0)", "pytest-mock (>=3.6.1)", "trustme", "uvloop (>=0.17)"] trio = ["trio (>=0.23)"] +[[package]] +name = "appnope" +version = "0.1.4" +description = "Disable App Nap on macOS >= 10.9" +optional = false +python-versions = ">=3.6" +files = [ + {file = "appnope-0.1.4-py2.py3-none-any.whl", hash = "sha256:502575ee11cd7a28c0205f379b525beefebab9d161b7c964670864014ed7213c"}, + {file = "appnope-0.1.4.tar.gz", hash = "sha256:1de3860566df9caf38f01f86f65e0e13e379af54f9e4bee1e66b48f2efffd1ee"}, +] + [[package]] name = "astroid" version = "3.1.0" @@ -1058,6 +1069,37 @@ toolz = ">=0.8.0" [package.extras] cython = ["cython"] +[[package]] +name = "debugpy" +version = "1.8.1" +description = "An implementation of the Debug Adapter Protocol for Python" +optional = false +python-versions = ">=3.8" +files = [ + {file = "debugpy-1.8.1-cp310-cp310-macosx_11_0_x86_64.whl", hash = "sha256:3bda0f1e943d386cc7a0e71bfa59f4137909e2ed947fb3946c506e113000f741"}, + {file = "debugpy-1.8.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:dda73bf69ea479c8577a0448f8c707691152e6c4de7f0c4dec5a4bc11dee516e"}, + {file = "debugpy-1.8.1-cp310-cp310-win32.whl", hash = "sha256:3a79c6f62adef994b2dbe9fc2cc9cc3864a23575b6e387339ab739873bea53d0"}, + {file = "debugpy-1.8.1-cp310-cp310-win_amd64.whl", hash = "sha256:7eb7bd2b56ea3bedb009616d9e2f64aab8fc7000d481faec3cd26c98a964bcdd"}, + {file = "debugpy-1.8.1-cp311-cp311-macosx_11_0_universal2.whl", hash = "sha256:016a9fcfc2c6b57f939673c874310d8581d51a0fe0858e7fac4e240c5eb743cb"}, + {file = "debugpy-1.8.1-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:fd97ed11a4c7f6d042d320ce03d83b20c3fb40da892f994bc041bbc415d7a099"}, + {file = "debugpy-1.8.1-cp311-cp311-win32.whl", hash = "sha256:0de56aba8249c28a300bdb0672a9b94785074eb82eb672db66c8144fff673146"}, + {file = "debugpy-1.8.1-cp311-cp311-win_amd64.whl", hash = "sha256:1a9fe0829c2b854757b4fd0a338d93bc17249a3bf69ecf765c61d4c522bb92a8"}, + {file = "debugpy-1.8.1-cp312-cp312-macosx_11_0_universal2.whl", hash = "sha256:3ebb70ba1a6524d19fa7bb122f44b74170c447d5746a503e36adc244a20ac539"}, + {file = "debugpy-1.8.1-cp312-cp312-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:a2e658a9630f27534e63922ebf655a6ab60c370f4d2fc5c02a5b19baf4410ace"}, + {file = "debugpy-1.8.1-cp312-cp312-win32.whl", hash = "sha256:caad2846e21188797a1f17fc09c31b84c7c3c23baf2516fed5b40b378515bbf0"}, + {file = "debugpy-1.8.1-cp312-cp312-win_amd64.whl", hash = "sha256:edcc9f58ec0fd121a25bc950d4578df47428d72e1a0d66c07403b04eb93bcf98"}, + {file = "debugpy-1.8.1-cp38-cp38-macosx_11_0_x86_64.whl", hash = "sha256:7a3afa222f6fd3d9dfecd52729bc2e12c93e22a7491405a0ecbf9e1d32d45b39"}, + {file = "debugpy-1.8.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d915a18f0597ef685e88bb35e5d7ab968964b7befefe1aaea1eb5b2640b586c7"}, + {file = "debugpy-1.8.1-cp38-cp38-win32.whl", hash = "sha256:92116039b5500633cc8d44ecc187abe2dfa9b90f7a82bbf81d079fcdd506bae9"}, + {file = "debugpy-1.8.1-cp38-cp38-win_amd64.whl", hash = "sha256:e38beb7992b5afd9d5244e96ad5fa9135e94993b0c551ceebf3fe1a5d9beb234"}, + {file = "debugpy-1.8.1-cp39-cp39-macosx_11_0_x86_64.whl", hash = "sha256:bfb20cb57486c8e4793d41996652e5a6a885b4d9175dd369045dad59eaacea42"}, + {file = "debugpy-1.8.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:efd3fdd3f67a7e576dd869c184c5dd71d9aaa36ded271939da352880c012e703"}, + {file = "debugpy-1.8.1-cp39-cp39-win32.whl", hash = "sha256:58911e8521ca0c785ac7a0539f1e77e0ce2df753f786188f382229278b4cdf23"}, + {file = "debugpy-1.8.1-cp39-cp39-win_amd64.whl", hash = "sha256:6df9aa9599eb05ca179fb0b810282255202a66835c6efb1d112d21ecb830ddd3"}, + {file = "debugpy-1.8.1-py2.py3-none-any.whl", hash = "sha256:28acbe2241222b87e255260c76741e1fbf04fdc3b6d094fcf57b6c6f75ce1242"}, + {file = "debugpy-1.8.1.zip", hash = "sha256:f696d6be15be87aef621917585f9bb94b1dc9e8aced570db1b8a6fc14e8f9b42"}, +] + [[package]] name = "decorator" version = "5.1.1" @@ -1307,53 +1349,53 @@ files = [ [[package]] name = "fonttools" -version = "4.52.4" +version = "4.53.0" description = "Tools to manipulate font files" optional = false python-versions = ">=3.8" files = [ - {file = "fonttools-4.52.4-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:fb8cd6559f0ae3a8f5e146f80ab2a90ad0325a759be8d48ee82758a0b89fa0aa"}, - {file = "fonttools-4.52.4-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:5ecb88318ff249bd2a715e7aec36774ce7ae3441128007ef72a39a60601f4a8f"}, - {file = "fonttools-4.52.4-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:b9a22cf1adaae7b2ba2ed7d8651a4193a4f348744925b4b740e6b38a94599c5b"}, - {file = "fonttools-4.52.4-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:8873d6edd1dae5c088dd3d61c9fd4dd80c827c486fa224d368233e7f33dc98af"}, - {file = "fonttools-4.52.4-cp310-cp310-musllinux_1_2_aarch64.whl", hash = "sha256:73ba38b98c012957940a04d9eb5439b42565ac892bba8cfc32e10d88e73921fe"}, - {file = "fonttools-4.52.4-cp310-cp310-musllinux_1_2_x86_64.whl", hash = "sha256:9725687db3c1cef13c0f40b380c3c15bea0113f4d0231b204d58edd5f2a53d90"}, - {file = "fonttools-4.52.4-cp310-cp310-win32.whl", hash = "sha256:9180775c9535389a665cae7c5282f8e07754beabf59b66aeba7f6bfeb32a3652"}, - {file = "fonttools-4.52.4-cp310-cp310-win_amd64.whl", hash = "sha256:46cc5d06ee05fd239c45d7935aaffd060ee773a88b97e901df50478247472643"}, - {file = "fonttools-4.52.4-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:d272c7e173c3085308345ccc7fb2ad6ce7f415d777791dd6ce4e8140e354d09c"}, - {file = "fonttools-4.52.4-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:21921e5855c399d10ddfc373538b425cabcf8b3258720b51450909e108896450"}, - {file = "fonttools-4.52.4-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:52f6001814ec5e0c961cabe89642f7e8d7e07892b565057aa526569b9ebb711c"}, - {file = "fonttools-4.52.4-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:4b0b9eb0f55dce9c7278ad4175f1cbaed23b799dce5ecc20e3213da241584140"}, - {file = "fonttools-4.52.4-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:70d87f2099006304d33438bdaa5101953b7e22e23a93b1c7b7ed0f32ff44b423"}, - {file = "fonttools-4.52.4-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:e176249292eccd89f81d39f514f2b5e8c75dfc9cef8653bdc3021d06697e9eff"}, - {file = "fonttools-4.52.4-cp311-cp311-win32.whl", hash = "sha256:bb7d206fa5ba6e082ba5d5e1b7107731029fc3a55c71c48de65121710d817986"}, - {file = "fonttools-4.52.4-cp311-cp311-win_amd64.whl", hash = "sha256:346d08ff92e577b2dc5a0c228487667d23fe2da35a8b9a8bba22c2b6ba8be21c"}, - {file = "fonttools-4.52.4-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:d2cc7906bc0afdd2689aaf88b910307333b1f936262d1d98f25dbf8a5eb2e829"}, - {file = "fonttools-4.52.4-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:00d9abf4b400f98fb895566eb298f60432b4b29048e3dc02807427b09a06604e"}, - {file = "fonttools-4.52.4-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:4b419207e53db1599b3d385afd4bca6692c219d53732890d0814a2593104d0e2"}, - {file = "fonttools-4.52.4-cp312-cp312-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:cf694159528022daa71b1777cb6ec9e0ebbdd29859f3e9c845826cafaef4ca29"}, - {file = "fonttools-4.52.4-cp312-cp312-musllinux_1_2_aarch64.whl", hash = "sha256:9a5d1b0475050056d2e3bc378014f2ea2230e8ae434eeac8dfb182aa8efaf642"}, - {file = "fonttools-4.52.4-cp312-cp312-musllinux_1_2_x86_64.whl", hash = "sha256:4c3ad89204c2d7f419436f1d6fde681b070c5e20b888beb57ccf92f640628cc9"}, - {file = "fonttools-4.52.4-cp312-cp312-win32.whl", hash = "sha256:1dc626de4b204d025d029e646bae8fdbf5acd9217158283a567f4b523fda3bae"}, - {file = "fonttools-4.52.4-cp312-cp312-win_amd64.whl", hash = "sha256:309b617942041073ffa96090d320b99d75648ed16e0c67fb1aa7788e06c834de"}, - {file = "fonttools-4.52.4-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:8b186cd6b8844f6cf04a7e0a174bc3649d3deddbfc10dc59846a4381f796d348"}, - {file = "fonttools-4.52.4-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:9ed23a03b7d9f0e29ca0679eafe5152aeccb0580312a3fc36f0662e178b4791b"}, - {file = "fonttools-4.52.4-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:89b53386214197bd5b3e3c753895bad691de84726ced3c222a59cde1dd12d57b"}, - {file = "fonttools-4.52.4-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:7467161f1eed557dbcec152d5ee95540200b1935709fa73307da16bc0b7ca361"}, - {file = "fonttools-4.52.4-cp38-cp38-musllinux_1_2_aarch64.whl", hash = "sha256:b4cba644e2515d685d4ee3ca2fbb5d53930a0e9ec2cf332ed704dc341b145878"}, - {file = "fonttools-4.52.4-cp38-cp38-musllinux_1_2_x86_64.whl", hash = "sha256:890e7a657574610330e42dd1e38d3b9e0a8cb0eff3da080f80995460a256d3dd"}, - {file = "fonttools-4.52.4-cp38-cp38-win32.whl", hash = "sha256:7dccf4666f716e5e0753f0fa28dad2f4431154c87747bc781c838b8a5dca990e"}, - {file = "fonttools-4.52.4-cp38-cp38-win_amd64.whl", hash = "sha256:a791f002d1b717268235cfae7e4957b7fd132e92e2c5400e521bf191f1b3a9a5"}, - {file = "fonttools-4.52.4-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:05e4291db6af66f466a203d9922e4c1d3e18ef16868f76f10b00e2c3b9814df2"}, - {file = "fonttools-4.52.4-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:a64e72d2c144630e017ac9c1c416ddf8ac43bef9a083bf81fe08c0695f0baa95"}, - {file = "fonttools-4.52.4-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:ebb183ed8b789cece0bd6363121913fb6da4034af89a2fa5408e42a1592889a8"}, - {file = "fonttools-4.52.4-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:a4daf2751a98c69d9620717826ed6c5743b662ef0ae7bb33dc6c205425e48eba"}, - {file = "fonttools-4.52.4-cp39-cp39-musllinux_1_2_aarch64.whl", hash = "sha256:15efb2ba4b8c2d012ee0bb7a850c2e4780c530cc83ec8e843b2a97f8b3a5fd4b"}, - {file = "fonttools-4.52.4-cp39-cp39-musllinux_1_2_x86_64.whl", hash = "sha256:35af630404223273f1d7acd4761f399131c62820366f53eac029337069f5826a"}, - {file = "fonttools-4.52.4-cp39-cp39-win32.whl", hash = "sha256:d0184aa88865339d96f7f452e8c5b621186ef7638744d78bf9b775d67e206819"}, - {file = "fonttools-4.52.4-cp39-cp39-win_amd64.whl", hash = "sha256:e03dae26084bb3632b4a77b1cd0419159d2226911aff6dc4c7e3058df68648c6"}, - {file = "fonttools-4.52.4-py3-none-any.whl", hash = "sha256:95e8a5975d08d0b624a14eec0f987e204ad81b480e24c5436af99170054434b8"}, - {file = "fonttools-4.52.4.tar.gz", hash = "sha256:859399b7adc8ac067be8e5c80ef4bb2faddff97e9b40896a9de75606a43d0469"}, + {file = "fonttools-4.53.0-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:52a6e0a7a0bf611c19bc8ec8f7592bdae79c8296c70eb05917fd831354699b20"}, + {file = "fonttools-4.53.0-cp310-cp310-macosx_11_0_arm64.whl", hash = "sha256:099634631b9dd271d4a835d2b2a9e042ccc94ecdf7e2dd9f7f34f7daf333358d"}, + {file = "fonttools-4.53.0-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:e40013572bfb843d6794a3ce076c29ef4efd15937ab833f520117f8eccc84fd6"}, + {file = "fonttools-4.53.0-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:715b41c3e231f7334cbe79dfc698213dcb7211520ec7a3bc2ba20c8515e8a3b5"}, + {file = "fonttools-4.53.0-cp310-cp310-musllinux_1_2_aarch64.whl", hash = "sha256:74ae2441731a05b44d5988d3ac2cf784d3ee0a535dbed257cbfff4be8bb49eb9"}, + {file = "fonttools-4.53.0-cp310-cp310-musllinux_1_2_x86_64.whl", hash = "sha256:95db0c6581a54b47c30860d013977b8a14febc206c8b5ff562f9fe32738a8aca"}, + {file = "fonttools-4.53.0-cp310-cp310-win32.whl", hash = "sha256:9cd7a6beec6495d1dffb1033d50a3f82dfece23e9eb3c20cd3c2444d27514068"}, + {file = "fonttools-4.53.0-cp310-cp310-win_amd64.whl", hash = "sha256:daaef7390e632283051e3cf3e16aff2b68b247e99aea916f64e578c0449c9c68"}, + {file = "fonttools-4.53.0-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:a209d2e624ba492df4f3bfad5996d1f76f03069c6133c60cd04f9a9e715595ec"}, + {file = "fonttools-4.53.0-cp311-cp311-macosx_11_0_arm64.whl", hash = "sha256:4f520d9ac5b938e6494f58a25c77564beca7d0199ecf726e1bd3d56872c59749"}, + {file = "fonttools-4.53.0-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:eceef49f457253000e6a2d0f7bd08ff4e9fe96ec4ffce2dbcb32e34d9c1b8161"}, + {file = "fonttools-4.53.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:fa1f3e34373aa16045484b4d9d352d4c6b5f9f77ac77a178252ccbc851e8b2ee"}, + {file = "fonttools-4.53.0-cp311-cp311-musllinux_1_2_aarch64.whl", hash = "sha256:28d072169fe8275fb1a0d35e3233f6df36a7e8474e56cb790a7258ad822b6fd6"}, + {file = "fonttools-4.53.0-cp311-cp311-musllinux_1_2_x86_64.whl", hash = "sha256:4a2a6ba400d386e904fd05db81f73bee0008af37799a7586deaa4aef8cd5971e"}, + {file = "fonttools-4.53.0-cp311-cp311-win32.whl", hash = "sha256:bb7273789f69b565d88e97e9e1da602b4ee7ba733caf35a6c2affd4334d4f005"}, + {file = "fonttools-4.53.0-cp311-cp311-win_amd64.whl", hash = "sha256:9fe9096a60113e1d755e9e6bda15ef7e03391ee0554d22829aa506cdf946f796"}, + {file = "fonttools-4.53.0-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:d8f191a17369bd53a5557a5ee4bab91d5330ca3aefcdf17fab9a497b0e7cff7a"}, + {file = "fonttools-4.53.0-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:93156dd7f90ae0a1b0e8871032a07ef3178f553f0c70c386025a808f3a63b1f4"}, + {file = "fonttools-4.53.0-cp312-cp312-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:bff98816cb144fb7b85e4b5ba3888a33b56ecef075b0e95b95bcd0a5fbf20f06"}, + {file = "fonttools-4.53.0-cp312-cp312-manylinux_2_5_x86_64.manylinux1_x86_64.manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:973d030180eca8255b1bce6ffc09ef38a05dcec0e8320cc9b7bcaa65346f341d"}, + {file = "fonttools-4.53.0-cp312-cp312-musllinux_1_2_aarch64.whl", hash = "sha256:c4ee5a24e281fbd8261c6ab29faa7fd9a87a12e8c0eed485b705236c65999109"}, + {file = "fonttools-4.53.0-cp312-cp312-musllinux_1_2_x86_64.whl", hash = "sha256:bd5bc124fae781a4422f61b98d1d7faa47985f663a64770b78f13d2c072410c2"}, + {file = "fonttools-4.53.0-cp312-cp312-win32.whl", hash = "sha256:a239afa1126b6a619130909c8404070e2b473dd2b7fc4aacacd2e763f8597fea"}, + {file = "fonttools-4.53.0-cp312-cp312-win_amd64.whl", hash = "sha256:45b4afb069039f0366a43a5d454bc54eea942bfb66b3fc3e9a2c07ef4d617380"}, + {file = "fonttools-4.53.0-cp38-cp38-macosx_10_9_universal2.whl", hash = "sha256:93bc9e5aaa06ff928d751dc6be889ff3e7d2aa393ab873bc7f6396a99f6fbb12"}, + {file = "fonttools-4.53.0-cp38-cp38-macosx_11_0_arm64.whl", hash = "sha256:2367d47816cc9783a28645bc1dac07f8ffc93e0f015e8c9fc674a5b76a6da6e4"}, + {file = "fonttools-4.53.0-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:907fa0b662dd8fc1d7c661b90782ce81afb510fc4b7aa6ae7304d6c094b27bce"}, + {file = "fonttools-4.53.0-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:3e0ad3c6ea4bd6a289d958a1eb922767233f00982cf0fe42b177657c86c80a8f"}, + {file = "fonttools-4.53.0-cp38-cp38-musllinux_1_2_aarch64.whl", hash = "sha256:73121a9b7ff93ada888aaee3985a88495489cc027894458cb1a736660bdfb206"}, + {file = "fonttools-4.53.0-cp38-cp38-musllinux_1_2_x86_64.whl", hash = "sha256:ee595d7ba9bba130b2bec555a40aafa60c26ce68ed0cf509983e0f12d88674fd"}, + {file = "fonttools-4.53.0-cp38-cp38-win32.whl", hash = "sha256:fca66d9ff2ac89b03f5aa17e0b21a97c21f3491c46b583bb131eb32c7bab33af"}, + {file = "fonttools-4.53.0-cp38-cp38-win_amd64.whl", hash = "sha256:31f0e3147375002aae30696dd1dc596636abbd22fca09d2e730ecde0baad1d6b"}, + {file = "fonttools-4.53.0-cp39-cp39-macosx_10_9_universal2.whl", hash = "sha256:7d6166192dcd925c78a91d599b48960e0a46fe565391c79fe6de481ac44d20ac"}, + {file = "fonttools-4.53.0-cp39-cp39-macosx_11_0_arm64.whl", hash = "sha256:ef50ec31649fbc3acf6afd261ed89d09eb909b97cc289d80476166df8438524d"}, + {file = "fonttools-4.53.0-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:7f193f060391a455920d61684a70017ef5284ccbe6023bb056e15e5ac3de11d1"}, + {file = "fonttools-4.53.0-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ba9f09ff17f947392a855e3455a846f9855f6cf6bec33e9a427d3c1d254c712f"}, + {file = "fonttools-4.53.0-cp39-cp39-musllinux_1_2_aarch64.whl", hash = "sha256:0c555e039d268445172b909b1b6bdcba42ada1cf4a60e367d68702e3f87e5f64"}, + {file = "fonttools-4.53.0-cp39-cp39-musllinux_1_2_x86_64.whl", hash = "sha256:5a4788036201c908079e89ae3f5399b33bf45b9ea4514913f4dbbe4fac08efe0"}, + {file = "fonttools-4.53.0-cp39-cp39-win32.whl", hash = "sha256:d1a24f51a3305362b94681120c508758a88f207fa0a681c16b5a4172e9e6c7a9"}, + {file = "fonttools-4.53.0-cp39-cp39-win_amd64.whl", hash = "sha256:1e677bfb2b4bd0e5e99e0f7283e65e47a9814b0486cb64a41adf9ef110e078f2"}, + {file = "fonttools-4.53.0-py3-none-any.whl", hash = "sha256:6b4f04b1fbc01a3569d63359f2227c89ab294550de277fd09d8fca6185669fa4"}, + {file = "fonttools-4.53.0.tar.gz", hash = "sha256:c93ed66d32de1559b6fc348838c7572d5c0ac1e4a258e76763a5caddd8944002"}, ] [package.extras] @@ -1806,6 +1848,39 @@ files = [ {file = "intel_openmp-2021.4.0-py2.py3-none-win_amd64.whl", hash = "sha256:eef4c8bcc8acefd7f5cd3b9384dbf73d59e2c99fc56545712ded913f43c4a94f"}, ] +[[package]] +name = "ipykernel" +version = "6.29.4" +description = "IPython Kernel for Jupyter" +optional = false +python-versions = ">=3.8" +files = [ + {file = "ipykernel-6.29.4-py3-none-any.whl", hash = "sha256:1181e653d95c6808039c509ef8e67c4126b3b3af7781496c7cbfb5ed938a27da"}, + {file = "ipykernel-6.29.4.tar.gz", hash = "sha256:3d44070060f9475ac2092b760123fadf105d2e2493c24848b6691a7c4f42af5c"}, +] + +[package.dependencies] +appnope = {version = "*", markers = "platform_system == \"Darwin\""} +comm = ">=0.1.1" +debugpy = ">=1.6.5" +ipython = ">=7.23.1" +jupyter-client = ">=6.1.12" +jupyter-core = ">=4.12,<5.0.dev0 || >=5.1.dev0" +matplotlib-inline = ">=0.1" +nest-asyncio = "*" +packaging = "*" +psutil = "*" +pyzmq = ">=24" +tornado = ">=6.1" +traitlets = ">=5.4.0" + +[package.extras] +cov = ["coverage[toml]", "curio", "matplotlib", "pytest-cov", "trio"] +docs = ["myst-parser", "pydata-sphinx-theme", "sphinx", "sphinx-autodoc-typehints", "sphinxcontrib-github-alt", "sphinxcontrib-spelling", "trio"] +pyqt5 = ["pyqt5"] +pyside6 = ["pyside6"] +test = ["flaky", "ipyparallel", "pre-commit", "pytest (>=7.0)", "pytest-asyncio (>=0.23.5)", "pytest-cov", "pytest-timeout"] + [[package]] name = "ipython" version = "8.18.1" @@ -2686,6 +2761,17 @@ nbformat = "*" sphinx = ">=1.8" traitlets = ">=5" +[[package]] +name = "nest-asyncio" +version = "1.6.0" +description = "Patch asyncio to allow nested event loops" +optional = false +python-versions = ">=3.5" +files = [ + {file = "nest_asyncio-1.6.0-py3-none-any.whl", hash = "sha256:87af6efd6b5e897c81050477ef65c62e2b2f35d51703cae01aff2905b1852e1c"}, + {file = "nest_asyncio-1.6.0.tar.gz", hash = "sha256:6f172d5449aca15afd6c646851f4e31e02c598d553a667e38cafa997cfec55fe"}, +] + [[package]] name = "networkx" version = "3.2.1" @@ -3878,6 +3964,7 @@ files = [ {file = "PyYAML-6.0.1-cp310-cp310-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:69b023b2b4daa7548bcfbd4aa3da05b3a74b772db9e23b982788168117739938"}, {file = "PyYAML-6.0.1-cp310-cp310-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:81e0b275a9ecc9c0c0c07b4b90ba548307583c125f54d5b6946cfee6360c733d"}, {file = "PyYAML-6.0.1-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:ba336e390cd8e4d1739f42dfe9bb83a3cc2e80f567d8805e11b46f4a943f5515"}, + {file = "PyYAML-6.0.1-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:326c013efe8048858a6d312ddd31d56e468118ad4cdeda36c719bf5bb6192290"}, {file = "PyYAML-6.0.1-cp310-cp310-win32.whl", hash = "sha256:bd4af7373a854424dabd882decdc5579653d7868b8fb26dc7d0e99f823aa5924"}, {file = "PyYAML-6.0.1-cp310-cp310-win_amd64.whl", hash = "sha256:fd1592b3fdf65fff2ad0004b5e363300ef59ced41c2e6b3a99d4089fa8c5435d"}, {file = "PyYAML-6.0.1-cp311-cp311-macosx_10_9_x86_64.whl", hash = "sha256:6965a7bc3cf88e5a1c3bd2e0b5c22f8d677dc88a455344035f03399034eb3007"}, @@ -3885,8 +3972,15 @@ files = [ {file = "PyYAML-6.0.1-cp311-cp311-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:42f8152b8dbc4fe7d96729ec2b99c7097d656dc1213a3229ca5383f973a5ed6d"}, {file = "PyYAML-6.0.1-cp311-cp311-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:062582fca9fabdd2c8b54a3ef1c978d786e0f6b3a1510e0ac93ef59e0ddae2bc"}, {file = "PyYAML-6.0.1-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:d2b04aac4d386b172d5b9692e2d2da8de7bfb6c387fa4f801fbf6fb2e6ba4673"}, + {file = "PyYAML-6.0.1-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:e7d73685e87afe9f3b36c799222440d6cf362062f78be1013661b00c5c6f678b"}, {file = "PyYAML-6.0.1-cp311-cp311-win32.whl", hash = "sha256:1635fd110e8d85d55237ab316b5b011de701ea0f29d07611174a1b42f1444741"}, {file = "PyYAML-6.0.1-cp311-cp311-win_amd64.whl", hash = "sha256:bf07ee2fef7014951eeb99f56f39c9bb4af143d8aa3c21b1677805985307da34"}, + {file = "PyYAML-6.0.1-cp312-cp312-macosx_10_9_x86_64.whl", hash = "sha256:855fb52b0dc35af121542a76b9a84f8d1cd886ea97c84703eaa6d88e37a2ad28"}, + {file = "PyYAML-6.0.1-cp312-cp312-macosx_11_0_arm64.whl", hash = "sha256:40df9b996c2b73138957fe23a16a4f0ba614f4c0efce1e9406a184b6d07fa3a9"}, + {file = "PyYAML-6.0.1-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:6c22bec3fbe2524cde73d7ada88f6566758a8f7227bfbf93a408a9d86bcc12a0"}, + {file = "PyYAML-6.0.1-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:8d4e9c88387b0f5c7d5f281e55304de64cf7f9c0021a3525bd3b1c542da3b0e4"}, + {file = "PyYAML-6.0.1-cp312-cp312-win32.whl", hash = "sha256:d483d2cdf104e7c9fa60c544d92981f12ad66a457afae824d146093b8c294c54"}, + {file = "PyYAML-6.0.1-cp312-cp312-win_amd64.whl", hash = "sha256:0d3304d8c0adc42be59c5f8a4d9e3d7379e6955ad754aa9d6ab7a398b59dd1df"}, {file = "PyYAML-6.0.1-cp36-cp36m-macosx_10_9_x86_64.whl", hash = "sha256:50550eb667afee136e9a77d6dc71ae76a44df8b3e51e41b77f6de2932bfe0f47"}, {file = "PyYAML-6.0.1-cp36-cp36m-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:1fe35611261b29bd1de0070f0b2f47cb6ff71fa6595c077e42bd0c419fa27b98"}, {file = "PyYAML-6.0.1-cp36-cp36m-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:704219a11b772aea0d8ecd7058d0082713c3562b4e271b849ad7dc4a5c90c13c"}, @@ -3903,6 +3997,7 @@ files = [ {file = "PyYAML-6.0.1-cp38-cp38-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:a0cd17c15d3bb3fa06978b4e8958dcdc6e0174ccea823003a106c7d4d7899ac5"}, {file = "PyYAML-6.0.1-cp38-cp38-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:28c119d996beec18c05208a8bd78cbe4007878c6dd15091efb73a30e90539696"}, {file = "PyYAML-6.0.1-cp38-cp38-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:7e07cbde391ba96ab58e532ff4803f79c4129397514e1413a7dc761ccd755735"}, + {file = "PyYAML-6.0.1-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:49a183be227561de579b4a36efbb21b3eab9651dd81b1858589f796549873dd6"}, {file = "PyYAML-6.0.1-cp38-cp38-win32.whl", hash = "sha256:184c5108a2aca3c5b3d3bf9395d50893a7ab82a38004c8f61c258d4428e80206"}, {file = "PyYAML-6.0.1-cp38-cp38-win_amd64.whl", hash = "sha256:1e2722cc9fbb45d9b87631ac70924c11d3a401b2d7f410cc0e3bbf249f2dca62"}, {file = "PyYAML-6.0.1-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:9eb6caa9a297fc2c2fb8862bc5370d0303ddba53ba97e71f08023b6cd73d16a8"}, @@ -3910,6 +4005,7 @@ files = [ {file = "PyYAML-6.0.1-cp39-cp39-manylinux_2_17_aarch64.manylinux2014_aarch64.whl", hash = "sha256:5773183b6446b2c99bb77e77595dd486303b4faab2b086e7b17bc6bef28865f6"}, {file = "PyYAML-6.0.1-cp39-cp39-manylinux_2_17_s390x.manylinux2014_s390x.whl", hash = "sha256:b786eecbdf8499b9ca1d697215862083bd6d2a99965554781d0d8d1ad31e13a0"}, {file = "PyYAML-6.0.1-cp39-cp39-manylinux_2_17_x86_64.manylinux2014_x86_64.whl", hash = "sha256:bc1bf2925a1ecd43da378f4db9e4f799775d6367bdb94671027b73b393a7c42c"}, + {file = "PyYAML-6.0.1-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:04ac92ad1925b2cff1db0cfebffb6ffc43457495c9b3c39d3fcae417d7125dc5"}, {file = "PyYAML-6.0.1-cp39-cp39-win32.whl", hash = "sha256:faca3bdcf85b2fc05d06ff3fbc1f83e1391b3e724afa3feba7d13eeab355484c"}, {file = "PyYAML-6.0.1-cp39-cp39-win_amd64.whl", hash = "sha256:510c9deebc5c0225e8c96813043e62b680ba2f9c50a08d3724c7f28a747d1486"}, {file = "PyYAML-6.0.1.tar.gz", hash = "sha256:bfdf460b1736c775f2ba9f6a92bca30bc2095067b8a9d77876d1fad6cc3b4a43"}, @@ -4306,30 +4402,50 @@ files = [ {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-macosx_10_9_universal2.whl", hash = "sha256:b42169467c42b692c19cf539c38d4602069d8c1505e97b86387fcf7afb766e1d"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-macosx_13_0_arm64.whl", hash = "sha256:07238db9cbdf8fc1e9de2489a4f68474e70dffcb32232db7c08fa61ca0c7c462"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl", hash = "sha256:fff3573c2db359f091e1589c3d7c5fc2f86f5bdb6f24252c2d8e539d4e45f412"}, + {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-manylinux_2_24_aarch64.whl", hash = "sha256:aa2267c6a303eb483de8d02db2871afb5c5fc15618d894300b88958f729ad74f"}, + {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-musllinux_1_1_i686.whl", hash = "sha256:840f0c7f194986a63d2c2465ca63af8ccbbc90ab1c6001b1978f05119b5e7334"}, + {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-musllinux_1_1_x86_64.whl", hash = "sha256:024cfe1fc7c7f4e1aff4a81e718109e13409767e4f871443cbff3dba3578203d"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-win32.whl", hash = "sha256:c69212f63169ec1cfc9bb44723bf2917cbbd8f6191a00ef3410f5a7fe300722d"}, {file = "ruamel.yaml.clib-0.2.8-cp310-cp310-win_amd64.whl", hash = "sha256:cabddb8d8ead485e255fe80429f833172b4cadf99274db39abc080e068cbcc31"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-macosx_10_9_universal2.whl", hash = "sha256:bef08cd86169d9eafb3ccb0a39edb11d8e25f3dae2b28f5c52fd997521133069"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-macosx_13_0_arm64.whl", hash = "sha256:b16420e621d26fdfa949a8b4b47ade8810c56002f5389970db4ddda51dbff248"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_24_x86_64.whl", hash = "sha256:25c515e350e5b739842fc3228d662413ef28f295791af5e5110b543cf0b57d9b"}, + {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-manylinux_2_24_aarch64.whl", hash = "sha256:1707814f0d9791df063f8c19bb51b0d1278b8e9a2353abbb676c2f685dee6afe"}, + {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-musllinux_1_1_i686.whl", hash = "sha256:46d378daaac94f454b3a0e3d8d78cafd78a026b1d71443f4966c696b48a6d899"}, + {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-musllinux_1_1_x86_64.whl", hash = "sha256:09b055c05697b38ecacb7ac50bdab2240bfca1a0c4872b0fd309bb07dc9aa3a9"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-win32.whl", hash = "sha256:53a300ed9cea38cf5a2a9b069058137c2ca1ce658a874b79baceb8f892f915a7"}, {file = "ruamel.yaml.clib-0.2.8-cp311-cp311-win_amd64.whl", hash = "sha256:c2a72e9109ea74e511e29032f3b670835f8a59bbdc9ce692c5b4ed91ccf1eedb"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-macosx_10_9_universal2.whl", hash = "sha256:ebc06178e8821efc9692ea7544aa5644217358490145629914d8020042c24aa1"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-macosx_13_0_arm64.whl", hash = "sha256:edaef1c1200c4b4cb914583150dcaa3bc30e592e907c01117c08b13a07255ec2"}, {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-manylinux_2_17_x86_64.manylinux2014_x86_64.manylinux_2_28_x86_64.whl", hash = "sha256:d176b57452ab5b7028ac47e7b3cf644bcfdc8cacfecf7e71759f7f51a59e5c92"}, + {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-manylinux_2_24_aarch64.whl", hash = "sha256:1dc67314e7e1086c9fdf2680b7b6c2be1c0d8e3a8279f2e993ca2a7545fecf62"}, + {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-musllinux_1_1_i686.whl", hash = "sha256:3213ece08ea033eb159ac52ae052a4899b56ecc124bb80020d9bbceeb50258e9"}, + {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-musllinux_1_1_x86_64.whl", hash = "sha256:aab7fd643f71d7946f2ee58cc88c9b7bfc97debd71dcc93e03e2d174628e7e2d"}, + {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-win32.whl", hash = "sha256:5c365d91c88390c8d0a8545df0b5857172824b1c604e867161e6b3d59a827eaa"}, + {file = "ruamel.yaml.clib-0.2.8-cp312-cp312-win_amd64.whl", hash = "sha256:1758ce7d8e1a29d23de54a16ae867abd370f01b5a69e1a3ba75223eaa3ca1a1b"}, {file = "ruamel.yaml.clib-0.2.8-cp36-cp36m-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:a5aa27bad2bb83670b71683aae140a1f52b0857a2deff56ad3f6c13a017a26ed"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-macosx_10_9_x86_64.whl", hash = "sha256:c58ecd827313af6864893e7af0a3bb85fd529f862b6adbefe14643947cfe2942"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-macosx_12_0_arm64.whl", hash = "sha256:f481f16baec5290e45aebdc2a5168ebc6d35189ae6fea7a58787613a25f6e875"}, + {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-manylinux_2_24_aarch64.whl", hash = "sha256:77159f5d5b5c14f7c34073862a6b7d34944075d9f93e681638f6d753606c6ce6"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:7f67a1ee819dc4562d444bbafb135832b0b909f81cc90f7aa00260968c9ca1b3"}, + {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-musllinux_1_1_i686.whl", hash = "sha256:4ecbf9c3e19f9562c7fdd462e8d18dd902a47ca046a2e64dba80699f0b6c09b7"}, + {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-musllinux_1_1_x86_64.whl", hash = "sha256:87ea5ff66d8064301a154b3933ae406b0863402a799b16e4a1d24d9fbbcbe0d3"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-win32.whl", hash = "sha256:75e1ed13e1f9de23c5607fe6bd1aeaae21e523b32d83bb33918245361e9cc51b"}, {file = "ruamel.yaml.clib-0.2.8-cp37-cp37m-win_amd64.whl", hash = "sha256:3f215c5daf6a9d7bbed4a0a4f760f3113b10e82ff4c5c44bec20a68c8014f675"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-macosx_10_9_x86_64.whl", hash = "sha256:1b617618914cb00bf5c34d4357c37aa15183fa229b24767259657746c9077615"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-macosx_12_0_arm64.whl", hash = "sha256:a6a9ffd280b71ad062eae53ac1659ad86a17f59a0fdc7699fd9be40525153337"}, + {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-manylinux_2_24_aarch64.whl", hash = "sha256:305889baa4043a09e5b76f8e2a51d4ffba44259f6b4c72dec8ca56207d9c6fe1"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:700e4ebb569e59e16a976857c8798aee258dceac7c7d6b50cab63e080058df91"}, + {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-musllinux_1_1_i686.whl", hash = "sha256:e2b4c44b60eadec492926a7270abb100ef9f72798e18743939bdbf037aab8c28"}, + {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-musllinux_1_1_x86_64.whl", hash = "sha256:e79e5db08739731b0ce4850bed599235d601701d5694c36570a99a0c5ca41a9d"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-win32.whl", hash = "sha256:955eae71ac26c1ab35924203fda6220f84dce57d6d7884f189743e2abe3a9fbe"}, {file = "ruamel.yaml.clib-0.2.8-cp38-cp38-win_amd64.whl", hash = "sha256:56f4252222c067b4ce51ae12cbac231bce32aee1d33fbfc9d17e5b8d6966c312"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-macosx_10_9_x86_64.whl", hash = "sha256:03d1162b6d1df1caa3a4bd27aa51ce17c9afc2046c31b0ad60a0a96ec22f8001"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-macosx_12_0_arm64.whl", hash = "sha256:bba64af9fa9cebe325a62fa398760f5c7206b215201b0ec825005f1b18b9bccf"}, + {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-manylinux_2_24_aarch64.whl", hash = "sha256:a1a45e0bb052edf6a1d3a93baef85319733a888363938e1fc9924cb00c8df24c"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-manylinux_2_5_x86_64.manylinux1_x86_64.whl", hash = "sha256:da09ad1c359a728e112d60116f626cc9f29730ff3e0e7db72b9a2dbc2e4beed5"}, + {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-musllinux_1_1_i686.whl", hash = "sha256:184565012b60405d93838167f425713180b949e9d8dd0bbc7b49f074407c5a8b"}, + {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-musllinux_1_1_x86_64.whl", hash = "sha256:a75879bacf2c987c003368cf14bed0ffe99e8e85acfa6c0bfffc21a090f16880"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-win32.whl", hash = "sha256:84b554931e932c46f94ab306913ad7e11bba988104c5cff26d90d03f68258cd5"}, {file = "ruamel.yaml.clib-0.2.8-cp39-cp39-win_amd64.whl", hash = "sha256:25ac8c08322002b06fa1d49d1646181f0b2c72f5cbc15a85e80b4c30a544bb15"}, {file = "ruamel.yaml.clib-0.2.8.tar.gz", hash = "sha256:beb2e0404003de9a4cab9753a8805a8fe9320ee6673136ed7f04255fe60bb512"}, @@ -4422,6 +4538,27 @@ dev = ["cython-lint (>=0.12.2)", "doit (>=0.36.0)", "mypy", "pycodestyle", "pyde doc = ["jupyterlite-pyodide-kernel", "jupyterlite-sphinx (>=0.12.0)", "jupytext", "matplotlib (>=3.5)", "myst-nb", "numpydoc", "pooch", "pydata-sphinx-theme (>=0.15.2)", "sphinx (>=5.0.0)", "sphinx-design (>=0.4.0)"] test = ["array-api-strict", "asv", "gmpy2", "hypothesis (>=6.30)", "mpmath", "pooch", "pytest", "pytest-cov", "pytest-timeout", "pytest-xdist", "scikit-umfpack", "threadpoolctl"] +[[package]] +name = "seaborn" +version = "0.13.2" +description = "Statistical data visualization" +optional = false +python-versions = ">=3.8" +files = [ + {file = "seaborn-0.13.2-py3-none-any.whl", hash = "sha256:636f8336facf092165e27924f223d3c62ca560b1f2bb5dff7ab7fad265361987"}, + {file = "seaborn-0.13.2.tar.gz", hash = "sha256:93e60a40988f4d65e9f4885df477e2fdaff6b73a9ded434c1ab356dd57eefff7"}, +] + +[package.dependencies] +matplotlib = ">=3.4,<3.6.1 || >3.6.1" +numpy = ">=1.20,<1.24.0 || >1.24.0" +pandas = ">=1.2" + +[package.extras] +dev = ["flake8", "flit", "mypy", "pandas-stubs", "pre-commit", "pytest", "pytest-cov", "pytest-xdist"] +docs = ["ipykernel", "nbconvert", "numpydoc", "pydata_sphinx_theme (==0.10.0rc2)", "pyyaml", "sphinx (<6.0.0)", "sphinx-copybutton", "sphinx-design", "sphinx-issues"] +stats = ["scipy (>=1.7)", "statsmodels (>=0.12)"] + [[package]] name = "setuptools" version = "69.5.1" @@ -5384,18 +5521,18 @@ files = [ [[package]] name = "zipp" -version = "3.19.0" +version = "3.19.1" description = "Backport of pathlib-compatible object wrapper for zip files" optional = false python-versions = ">=3.8" files = [ - {file = "zipp-3.19.0-py3-none-any.whl", hash = "sha256:96dc6ad62f1441bcaccef23b274ec471518daf4fbbc580341204936a5a3dddec"}, - {file = "zipp-3.19.0.tar.gz", hash = "sha256:952df858fb3164426c976d9338d3961e8e8b3758e2e059e0f754b8c4262625ee"}, + {file = "zipp-3.19.1-py3-none-any.whl", hash = "sha256:2828e64edb5386ea6a52e7ba7cdb17bb30a73a858f5eb6eb93d8d36f5ea26091"}, + {file = "zipp-3.19.1.tar.gz", hash = "sha256:35427f6d5594f4acf82d25541438348c26736fa9b3afa2754bcd63cdb99d8e8f"}, ] [package.extras] -docs = ["furo", "jaraco.packaging (>=9.3)", "jaraco.tidelift (>=1.4)", "rst.linker (>=1.9)", "sphinx (>=3.5)", "sphinx-lint"] -testing = ["big-O", "jaraco.functools", "jaraco.itertools", "jaraco.test", "more-itertools", "pytest (>=6,!=8.1.*)", "pytest-checkdocs (>=2.4)", "pytest-cov", "pytest-enabler (>=2.2)", "pytest-ignore-flaky", "pytest-mypy", "pytest-ruff (>=0.2.1)"] +doc = ["furo", "jaraco.packaging (>=9.3)", "jaraco.tidelift (>=1.4)", "rst.linker (>=1.9)", "sphinx (>=3.5)", "sphinx-lint"] +test = ["big-O", "jaraco.functools", "jaraco.itertools", "jaraco.test", "more-itertools", "pytest (>=6,!=8.1.*)", "pytest-checkdocs (>=2.4)", "pytest-cov", "pytest-enabler (>=2.2)", "pytest-ignore-flaky", "pytest-mypy", "pytest-ruff (>=0.2.1)"] [extras] qinfo = [] @@ -5405,4 +5542,4 @@ torch = ["torch"] [metadata] lock-version = "2.0" python-versions = ">=3.9,<3.13" -content-hash = "5a79a55221a1417404521c98624455fac4eea5dcdac0b805a8f15c768bd29b65" +content-hash = "8e08af9af1dbb26d97fcb75be806be99ea4cc2619a6c986c48f627d3b7d9541a" diff --git a/pyproject.toml b/pyproject.toml index 2e76b0ea33..40a038f9ac 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -57,6 +57,8 @@ sphinx-copybutton = "^0.5.2" nbsphinx = "^0.8.12" ipython = "^8.10.0" qulacs = "^0.6.3" +seaborn = "^0.13.2" +ipykernel = "^6.29.4" [tool.poetry.group.tests] optional = true diff --git a/src/qibo/models/dbi/double_bracket.py b/src/qibo/models/dbi/double_bracket.py index 0c61692650..db3806c566 100644 --- a/src/qibo/models/dbi/double_bracket.py +++ b/src/qibo/models/dbi/double_bracket.py @@ -1,13 +1,18 @@ import math from copy import deepcopy from enum import Enum, auto -from functools import partial from typing import Optional -import hyperopt import numpy as np from qibo.hamiltonians import Hamiltonian +from qibo.models.dbi.utils import * +from qibo.models.dbi.utils_scheduling import ( + grid_search_step, + hyperopt_step, + polynomial_step, + simulated_annealing_step, +) class DoubleBracketGeneratorType(Enum): @@ -19,7 +24,36 @@ class DoubleBracketGeneratorType(Enum): """Use single commutator.""" group_commutator = auto() """Use group commutator approximation""" - # TODO: add double commutator (does it converge?) + group_commutator_third_order = auto() + """Implements: $e^{\frac{\\sqrt{5}-1}{2}isH}e^{\frac{\\sqrt{5}-1}{2}isD}e^{-isH}e^{isD}e^{\frac{3-\\sqrt{5}}{2}isH}e^{isD} + \approx e^{-s^2[H,D]} + O(s^4)$ + which is equation (8) in https://arxiv.org/abs/2111.12177] + s must be taken as $\\sqrt{s}$ to approximate the flow using the commutator + """ + + +class DoubleBracketCostFunction(str, Enum): + """Define the DBI cost function.""" + + off_diagonal_norm = "off_diagonal_norm" + """Use off-diagonal norm as cost function.""" + least_squares = "least_squares" + """Use least squares as cost function.""" + energy_fluctuation = "energy_fluctuation" + """Use energy fluctuation as cost function.""" + + +class DoubleBracketScheduling(Enum): + """Define the DBI scheduling strategies.""" + + hyperopt = hyperopt_step + """Use hyperopt package.""" + grid_search = grid_search_step + """Use greedy grid search.""" + polynomial_approximation = polynomial_step + """Use polynomial expansion (analytical) of the loss function.""" + simulated_annealing = simulated_annealing_step + """Use simulated annealing algorithm""" class DoubleBracketScheduling(Enum): @@ -61,16 +95,29 @@ def __init__( self, hamiltonian: Hamiltonian, mode: DoubleBracketGeneratorType = DoubleBracketGeneratorType.canonical, - scheduling: DoubleBracketScheduling = DoubleBracketScheduling.use_grid_search, + scheduling: DoubleBracketScheduling = DoubleBracketScheduling.grid_search, + cost: DoubleBracketCostFunction = DoubleBracketCostFunction.off_diagonal_norm, + ref_state: np.array = None, ): self.h = hamiltonian self.h0 = deepcopy(self.h) self.mode = mode self.scheduling = scheduling + self.cost = cost + self.ref_state = ref_state + """ + Args: + hamiltonian (Hamiltonian): Starting Hamiltonian; + mode (DoubleBracketGeneratorType): type of generator of the evolution. + scheduling (DoubleBracketScheduling): type of scheduling strategy. + cost (DoubleBracketCost): type of cost function. + ref_state (np.array): reference state for computing the energy fluctuation. + """ def __call__( self, step: float, mode: DoubleBracketGeneratorType = None, d: np.array = None ): + """Performs one double bracket rotation.""" if mode is None: mode = self.mode @@ -83,8 +130,8 @@ def __call__( if d is None: d = self.diagonal_h_matrix operator = self.backend.calculate_matrix_exp( - -1.0j * step, - self.commutator(d, self.h.matrix), + 1.0j * step, + self.commutator(self.backend.cast(d), self.h.matrix), ) elif mode is DoubleBracketGeneratorType.group_commutator: if d is None: @@ -95,7 +142,22 @@ def __call__( @ self.h.exp(step) @ self.backend.calculate_matrix_exp(step, d) ) - return operator + elif mode is DoubleBracketGeneratorType.group_commutator_third_order: + if d is None: + d = self.diagonal_h_matrix + operator = ( + self.h.exp(-step * (np.sqrt(5) - 1) / 2) + @ self.backend.calculate_matrix_exp(-step * (np.sqrt(5) - 1) / 2, d) + @ self.h.exp(step) + @ self.backend.calculate_matrix_exp(step * (np.sqrt(5) + 1) / 2, d) + @ self.h.exp(-step * (3 - np.sqrt(5)) / 2) + @ self.backend.calculate_matrix_exp(-step, d) + ) + operator_dagger = self.backend.cast( + np.array(np.matrix(self.backend.to_numpy(operator)).getH()) + ) + + self.h.matrix = operator @ self.h.matrix @ operator_dagger @staticmethod def commutator(a, b): @@ -109,11 +171,12 @@ def diagonal_h_matrix(self): @property def off_diag_h(self): + """Off-diagonal H matrix.""" return self.h.matrix - self.diagonal_h_matrix @property def off_diagonal_norm(self): - r"""Hilbert Schmidt norm of off-diagonal part of H matrix, namely :math:`\\text{Tr}(\\sqrt{A^{\\dagger} A})`.""" + """Hilbert Schmidt norm of off-diagonal part of H matrix, namely :math:`\\text{Tr}(\\sqrt{A^{\\dagger} A})`.""" off_diag_h_dag = self.backend.cast( np.matrix(self.backend.to_numpy(self.off_diag_h)).getH() ) @@ -126,79 +189,40 @@ def backend(self): """Get Hamiltonian's backend.""" return self.h0.backend - def grid_search_step( - self, - step_min: float = 1e-5, - step_max: float = 1, - num_evals: int = 100, - space: Optional[np.array] = None, - d: Optional[np.array] = None, - ): - """ - Greedy optimization of the iteration step. - - Args: - step_min: lower bound of the search grid; - step_max: upper bound of the search grid; - mnum_evals: number of iterations between step_min and step_max; - d: diagonal operator for generating double-bracket iterations. - - Returns: - (float): optimized best iteration step (minimizing off-diagonal norm). - """ - if space is None: - space = np.linspace(step_min, step_max, num_evals) - - if d is None: - d = self.diagonal_h_matrix - - loss_list = [self.loss(step, d=d) for step in space] - idx_max_loss = loss_list.index(min(loss_list)) - return space[idx_max_loss] + @property + def nqubits(self): + """Number of qubits.""" + return self.h.nqubits + + def least_squares(self, d: np.array): + """Least squares cost function.""" + d = self.backend.to_numpy(d) + return np.real( + 0.5 * np.linalg.norm(d) ** 2 + - np.trace(self.backend.to_numpy(self.h.matrix) @ d) + ) - def hyperopt_step( + def choose_step( self, - step_min: float = 1e-5, - step_max: float = 1, - max_evals: int = 1000, - space: callable = None, - optimizer: callable = None, - look_ahead: int = 1, - verbose: bool = False, d: Optional[np.array] = None, + scheduling: Optional[DoubleBracketScheduling] = None, + **kwargs, ): - """ - Optimize iteration step using hyperopt. - - Args: - step_min: lower bound of the search grid; - step_max: upper bound of the search grid; - max_evals: maximum number of iterations done by the hyperoptimizer; - space: see hyperopt.hp possibilities; - optimizer: see hyperopt algorithms; - look_ahead: number of iteration steps to compute the loss function; - verbose: level of verbosity; - d: diagonal operator for generating double-bracket iterations. - - Returns: - (float): optimized best iteration step (minimizing off-diagonal norm). - """ - if space is None: - space = hyperopt.hp.uniform - if optimizer is None: - optimizer = hyperopt.tpe - if d is None: - d = self.diagonal_h_matrix - - space = space("step", step_min, step_max) - best = hyperopt.fmin( - fn=partial(self.loss, d=d, look_ahead=look_ahead), - space=space, - algo=optimizer.suggest, - max_evals=max_evals, - verbose=verbose, - ) - return best["step"] + """Calculate the optimal step using respective the `scheduling` methods.""" + if scheduling is None: + scheduling = self.scheduling + step = scheduling(self, d=d, **kwargs) + # TODO: write test for this case + if ( + step is None + and scheduling is DoubleBracketScheduling.polynomial_approximation + ): # pragma: no cover + kwargs["n"] = kwargs.get("n", 3) + kwargs["n"] += 1 + # if n==n_max, return None + step = scheduling(self, d=d, **kwargs) + # if for a given polynomial order n, no solution is found, we increase the order of the polynomial by 1 + return step def polynomial_step( self, @@ -292,9 +316,9 @@ def loss(self, step: float, d: np.array = None, look_ahead: int = 1): Compute loss function distance between `look_ahead` steps. Args: - step: iteration step. - d: diagonal operator, use canonical by default. - look_ahead: number of iteration steps to compute the loss function; + step (float): iteration step. + d (np.array): diagonal operator, use canonical by default. + look_ahead (int): number of iteration steps to compute the loss function; """ # copy initial hamiltonian h_copy = deepcopy(self.h) @@ -302,8 +326,13 @@ def loss(self, step: float, d: np.array = None, look_ahead: int = 1): for _ in range(look_ahead): self.__call__(mode=self.mode, step=step, d=d) - # off_diagonal_norm's value after the steps - loss = self.off_diagonal_norm + # loss values depending on the cost function + if self.cost is DoubleBracketCostFunction.off_diagonal_norm: + loss = self.off_diagonal_norm + elif self.cost is DoubleBracketCostFunction.least_squares: + loss = self.least_squares(d) + elif self.cost == DoubleBracketCostFunction.energy_fluctuation: + loss = self.energy_fluctuation(self.ref_state) # set back the initial configuration self.h = h_copy @@ -312,10 +341,10 @@ def loss(self, step: float, d: np.array = None, look_ahead: int = 1): def energy_fluctuation(self, state): """ - Evaluate energy fluctuation + Evaluate energy fluctuation. .. math:: - \\Xi_{k}(\\mu) = \\sqrt{\\langle\\mu|\\hat{H}^2|\\mu\\rangle - \\langle\\mu|\\hat{H}|\\mu\\rangle^2} \\, + \\Xi(\\mu) = \\sqrt{\\langle\\mu|\\hat{H}^2|\\mu\\rangle - \\langle\\mu|\\hat{H}|\\mu\\rangle^2} \\, for a given state :math:`|\\mu\\rangle`. @@ -323,3 +352,32 @@ def energy_fluctuation(self, state): state (np.ndarray): quantum state to be used to compute the energy fluctuation with H. """ return self.h.energy_fluctuation(state) + + def sigma(self, h: np.array): + """Returns the off-diagonal restriction of matrix `h`.""" + return self.backend.cast(h) - self.backend.cast( + np.diag(np.diag(self.backend.to_numpy(h))) + ) + + def generate_gamma_list(self, n: int, d: np.array): + r"""Computes the n-nested Gamma functions, where $\Gamma_k=[W,...,[W,[W,H]]...]$, where we take k nested commutators with $W = [D, H]$""" + W = self.commutator(self.backend.cast(d), self.sigma(self.h.matrix)) + gamma_list = [self.h.matrix] + for _ in range(n - 1): + gamma_list.append(self.commutator(W, gamma_list[-1])) + return gamma_list + + def cost_expansion(self, d, n): + d = self.backend.cast(d) + + if self.cost is DoubleBracketCostFunction.off_diagonal_norm: + coef = off_diagonal_norm_polynomial_expansion_coef(self, d, n) + elif self.cost is DoubleBracketCostFunction.least_squares: + coef = least_squares_polynomial_expansion_coef(self, d, n) + elif self.cost is DoubleBracketCostFunction.energy_fluctuation: + coef = energy_fluctuation_polynomial_expansion_coef( + self, d, n, self.ref_state + ) + else: # pragma: no cover + raise ValueError(f"Cost function {self.cost} not recognized.") + return coef diff --git a/src/qibo/models/dbi/utils.py b/src/qibo/models/dbi/utils.py index 79bb41db7c..10367d71f9 100644 --- a/src/qibo/models/dbi/utils.py +++ b/src/qibo/models/dbi/utils.py @@ -1,17 +1,12 @@ -from copy import deepcopy -from itertools import product -from typing import Optional +import math +from enum import Enum, auto +from itertools import combinations, product import numpy as np -from hyperopt import hp, tpe from qibo import symbols from qibo.backends import _check_backend from qibo.hamiltonians import SymbolicHamiltonian -from qibo.models.dbi.double_bracket import ( - DoubleBracketGeneratorType, - DoubleBracketIteration, -) def generate_Z_operators(nqubits: int, backend=None): @@ -71,96 +66,200 @@ def str_to_symbolic(name: str): return tensor_op -def select_best_dbr_generator( - dbi_object: DoubleBracketIteration, - d_list: list, - step: Optional[float] = None, - step_min: float = 1e-5, - step_max: float = 1, - max_evals: int = 200, - compare_canonical: bool = True, +def cs_angle_sgn(dbi_object, d): + """Calculates the sign of Cauchy-Schwarz Angle :math:`\\langle W(Z), W({\\rm canonical}) \\rangle_{\\rm HS}`.""" + d = dbi_object.backend.cast(d) + norm = np.trace( + np.dot( + np.conjugate( + dbi_object.commutator(dbi_object.diagonal_h_matrix, dbi_object.h.matrix) + ).T, + dbi_object.commutator(d, dbi_object.h.matrix), + ) + ) + return np.real(np.sign(norm)) + + +def decompose_into_pauli_basis(h_matrix: np.array, pauli_operators: list): + """finds the decomposition of hamiltonian `h_matrix` into Pauli-Z operators""" + nqubits = int(np.log2(h_matrix.shape[0])) + + decomposition = [] + for Z_i in pauli_operators: + expect = np.trace(h_matrix @ Z_i) / 2**nqubits + decomposition.append(expect) + return decomposition + + +def generate_pauli_index(nqubits, order): + """ + Generate all possible combinations of qubits for a given order of Pauli operators. + """ + if order == 1: + return list(range(nqubits)) + else: + indices = list(range(nqubits)) + return indices + [ + comb for i in range(2, order + 1) for comb in combinations(indices, i) + ] + + +def generate_pauli_operator_dict( + nqubits: int, + parameterization_order: int = 1, + symbols_pauli=symbols.Z, + backend=None, ): - """Selects the best double bracket rotation generator from a list and execute the rotation. + """Generates a dictionary containing Pauli `symbols_pauli` operators of locality `parameterization_order` for `nqubits` qubits. Args: - dbi_object (`DoubleBracketIteration`): the target DoubleBracketIteration object. - d_list (list): list of diagonal operators (np.array) to run from. - step (float): fixed iteration duration. - Defaults to ``None``, uses hyperopt. - step_min (float): minimally allowed iteration duration. - step_max (float): maximally allowed iteration duration. - max_evals (int): maximally allowed number of evaluation in hyperopt. - compare_canonical (bool): if `True`, the optimal diagonal operator chosen from "d_list" is compared with the canonical bracket. + nqubits (int): number of qubits in the system. + parameterization_order (int, optional): the locality of the operators generated. Defaults to 1. + symbols_pauli (qibo.symbols, optional): the symbol of the intended Pauli operator. Defaults to symbols.Z. Returns: - The updated dbi_object, index of the optimal diagonal operator, respective step duration, and evolution direction. + pauli_operator_dict (dictionary): dictionary with structure {"operator_name": operator} + + Example: + pauli_operator_dict = generate_pauli_operator_dict) """ - norms_off_diagonal_restriction = [ - dbi_object.off_diagonal_norm for _ in range(len(d_list)) + backend = _check_backend(backend) + pauli_index = generate_pauli_index(nqubits, order=parameterization_order) + pauli_operators = [ + generate_pauli_operators(nqubits, symbols_pauli, index, backend=backend) + for index in pauli_index ] - optimal_steps, flip_list = [], [] - for i, d in enumerate(d_list): - # prescribed step durations - dbi_eval = deepcopy(dbi_object) - flip_list.append(cs_angle_sgn(dbi_eval, d)) - if flip_list[i] != 0: - if step is None: - step_best = dbi_eval.hyperopt_step( - d=flip_list[i] * d, - step_min=step_min, - step_max=step_max, - space=hp.uniform, - optimizer=tpe, - max_evals=max_evals, - ) - else: - step_best = step - dbi_eval(step=step_best, d=flip_list[i] * d) - optimal_steps.append(step_best) - norms_off_diagonal_restriction[i] = dbi_eval.off_diagonal_norm - # canonical - if compare_canonical is True: - flip_list.append(1) - dbi_eval = deepcopy(dbi_object) - dbi_eval.mode = DoubleBracketGeneratorType.canonical - if step is None: - step_best = dbi_eval.hyperopt_step( - step_min=step_min, - step_max=step_max, - space=hp.uniform, - optimizer=tpe, - max_evals=max_evals, - ) - else: - step_best = step - dbi_object(step=step) - optimal_steps.append(step) - norms_off_diagonal_restriction.append(dbi_object.off_diagonal_norm) - # find best d - idx_max_loss = norms_off_diagonal_restriction.index( - min(norms_off_diagonal_restriction) - ) - flip = flip_list[idx_max_loss] - step_optimal = optimal_steps[idx_max_loss] - dbi_eval = deepcopy(dbi_object) - if idx_max_loss == len(d_list) and compare_canonical is True: - # canonical - dbi_eval(step=step_optimal, mode=DoubleBracketGeneratorType.canonical) + return {index: operator for index, operator in zip(pauli_index, pauli_operators)} + +def generate_pauli_operators(nqubits, symbols_pauli, positions, backend=None): + # generate matrix of an nqubit-pauli operator with `symbols_pauli` at `positions` + if isinstance(positions, int): + return SymbolicHamiltonian( + symbols_pauli(positions), + nqubits=nqubits, + backend=backend, + ).dense.matrix else: - d_optimal = flip * d_list[idx_max_loss] - dbi_eval(step=step_optimal, d=d_optimal) - return dbi_eval, idx_max_loss, step_optimal, flip + terms = [symbols_pauli(pos) for pos in positions] + return SymbolicHamiltonian( + math.prod(terms), nqubits=nqubits, backend=backend + ).dense.matrix -def cs_angle_sgn(dbi_object, d): - """Calculates the sign of Cauchy-Schwarz Angle :math:`\\langle W(Z), W({\\rm canonical}) \\rangle_{\\rm HS}`.""" - norm = np.trace( - np.dot( - np.conjugate( - dbi_object.commutator(dbi_object.diagonal_h_matrix, dbi_object.h.matrix) - ).T, - dbi_object.commutator(d, dbi_object.h.matrix), +class ParameterizationTypes(Enum): + """Define types of parameterization for diagonal operator.""" + + pauli = auto() + """Uses Pauli-Z operators (magnetic field).""" + computational = auto() + """Uses computational basis.""" + + +def params_to_diagonal_operator( + params: np.array, + nqubits: int, + parameterization: ParameterizationTypes = ParameterizationTypes.pauli, + pauli_parameterization_order: int = 1, + normalize: bool = False, + pauli_operator_dict: dict = None, + backend=None, +): + r"""Creates the $D$ operator for the double-bracket iteration ansatz depending on the parameterization type.""" + backend = _check_backend(backend) + if parameterization is ParameterizationTypes.pauli: + # raise error if dimension mismatch + d = sum( + [ + backend.to_numpy(params[i]) + * backend.to_numpy(list(pauli_operator_dict.values())[i]) + for i in range(nqubits) + ] + ) + elif parameterization is ParameterizationTypes.computational: + d = np.zeros((len(params), len(params))) + for i in range(len(params)): + d[i, i] = backend.to_numpy(params[i]) + + # TODO: write proper tests for normalize=True + if normalize: # pragma: no cover + d = d / np.linalg.norm(d) + return d + + +def off_diagonal_norm_polynomial_expansion_coef(dbi_object, d, n): + # generate Gamma's where $\Gamma_{k+1}=[W, \Gamma_{k}], $\Gamma_0=H + W = dbi_object.commutator( + dbi_object.backend.cast(d), dbi_object.sigma(dbi_object.h.matrix) + ) + gamma_list = dbi_object.generate_gamma_list(n + 2, d) + sigma_gamma_list = list(map(dbi_object.sigma, gamma_list)) + gamma_list_np = list(map(dbi_object.backend.to_numpy, sigma_gamma_list)) + exp_list = np.array([1 / math.factorial(k) for k in range(n + 1)]) + # coefficients for rotation with [W,H] and H + c1 = exp_list.reshape((-1, 1, 1)) * gamma_list_np[1:] + c2 = exp_list.reshape((-1, 1, 1)) * gamma_list_np[:-1] + # product coefficient + trace_coefficients = [0] * (2 * n + 1) + for k in range(n + 1): + for j in range(n + 1): + power = k + j + product_matrix = c1[k] @ c2[j] + trace_coefficients[power] += 2 * np.trace(product_matrix) + # coefficients from high to low (n:0) + coef = list(reversed(trace_coefficients[: n + 1])) + return coef + + +def least_squares_polynomial_expansion_coef(dbi_object, d, n: int = 3): + """Return the Taylor expansion coefficients of least square cost of `dbi_object.h` and diagonal operator `d` with respect to double bracket rotation duration `s`.""" + # generate Gamma's where $\Gamma_{k+1}=[W, \Gamma_{k}], $\Gamma_0=H + Gamma_list = dbi_object.generate_gamma_list(n + 1, d) + exp_list = np.array([1 / math.factorial(k) for k in range(n + 1)]) + # coefficients + coef = np.empty(n) + for i in range(n): + coef[i] = np.real( + exp_list[i] * np.trace(dbi_object.backend.cast(d) @ Gamma_list[i + 1]) ) + coef = list(reversed(coef)) + return coef + + +def energy_fluctuation_polynomial_expansion_coef( + dbi_object, d: np.array, n: int = 3, state=0 +): + """Return the Taylor expansion coefficients of energy fluctuation of `dbi_object` with respect to double bracket rotation duration `s`.""" + # generate Gamma's where $\Gamma_{k+1}=[W, \Gamma_{k}], $\Gamma_0=H + Gamma_list = dbi_object.generate_gamma_list(n + 1, d) + # coefficients + coef = np.empty(3) + state_cast = dbi_object.backend.cast(state) + state_dag = dbi_object.backend.cast(state.conj().T) + + def variance(a): + """Calculates the variance of a matrix A with respect to a state: + Var($A$) = $\\langle\\mu|A^2|\\mu\rangle-\\langle\\mu|A|\\mu\rangle^2$""" + b = a @ a + return state_dag @ b @ state_cast - (state_dag @ a @ state_cast) ** 2 + + def covariance(a, b): + """This is a generalization of the notion of covariance, needed for the polynomial expansion of the energy fluctuation, + applied to two operators A and B with respect to a state: + Cov($A,B$) = $\\langle\\mu|AB|\\mu\rangle-\\langle\\mu|A|\\mu\rangle\\langle\\mu|B|\\mu\rangle$ + """ + + c = a @ b + b @ a + return ( + state_dag @ c @ state_cast + - 2 * state_dag @ a @ state_cast * state_dag @ b @ state_cast + ) + + coef[0] = np.real(2 * covariance(Gamma_list[0], Gamma_list[1])) + coef[1] = np.real(2 * variance(Gamma_list[1])) + coef[2] = np.real( + covariance(Gamma_list[0], Gamma_list[3]) + + 3 * covariance(Gamma_list[1], Gamma_list[2]) ) - return np.sign(norm) + coef = list(reversed(coef)) + return coef diff --git a/src/qibo/models/dbi/utils_dbr_strategies.py b/src/qibo/models/dbi/utils_dbr_strategies.py new file mode 100644 index 0000000000..5aae761fde --- /dev/null +++ b/src/qibo/models/dbi/utils_dbr_strategies.py @@ -0,0 +1,299 @@ +import hyperopt + +from qibo.backends import _check_backend +from qibo.models.dbi.double_bracket import * +from qibo.models.dbi.utils import * + + +def select_best_dbr_generator( + dbi_object: DoubleBracketIteration, + d_list: list, + step: Optional[float] = None, + compare_canonical: bool = True, + scheduling: DoubleBracketScheduling = None, + **kwargs, +): + """Selects the best double bracket rotation generator from a list and execute the rotation. + + Args: + dbi_object (`DoubleBracketIteration`): the target DoubleBracketIteration object. + d_list (list): list of diagonal operators (np.array) to select from. + step (float): fixed iteration duration. + Defaults to ``None``, optimize with `scheduling` method and `choose_step` function. + compare_canonical (boolean): if `True`, the diagonalization effect with operators from `d_list` is compared with the canonical bracket. + scheduling (`DoubleBracketScheduling`): scheduling method for finding the optimal step. + + Returns: + The updated dbi_object (`DoubleBracketIteration`), index of the optimal diagonal operator (int), respective step duration (float), and sign (int). + + Example: + from qibo.hamiltonians import Hamiltonian + from qibo.models.dbi.double_bracket import * + from qibo.models.dbi.utils_dbr_strategies import select_best_dbr_generator + from qibo.quantum_info import random_hermitian + + nqubits = 3 + NSTEPS = 3 + h0 = random_hermitian(2**nqubits) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0), + mode=DoubleBracketGeneratorType.single_commutator, + ) + initial_off_diagonal_norm = dbi.off_diagonal_norm + generate_local_Z = generate_Z_operators(nqubits) + Z_ops = list(generate_local_Z.values()) + for _ in range(NSTEPS): + dbi, idx, step, flip_sign = select_best_dbr_generator( + dbi, Z_ops, compare_canonical=True + ) + """ + if scheduling is None: + scheduling = dbi_object.scheduling + + if compare_canonical: + norms_off_diagonal_restriction = [dbi_object.off_diagonal_norm] * ( + len(d_list) + 1 + ) + optimal_steps = np.zeros(len(d_list) + 1) + flip_list = np.ones(len(d_list) + 1) + else: + norms_off_diagonal_restriction = [dbi_object.off_diagonal_norm] * (len(d_list)) + optimal_steps = np.zeros(len(d_list)) + flip_list = np.ones(len(d_list)) + + for i, d in enumerate(d_list): + # prescribed step durations + dbi_eval = deepcopy(dbi_object) + d = dbi_eval.backend.cast(d) + flip_list[i] = cs_angle_sgn(dbi_eval, d) + if flip_list[i] != 0: + if step is None: + step_best = dbi_eval.choose_step( + d=flip_list[i] * d, scheduling=scheduling, **kwargs + ) + else: + step_best = step + dbi_eval(step=step_best, d=flip_list[i] * d) + optimal_steps[i] = step_best + norms_off_diagonal_restriction[i] = dbi_eval.off_diagonal_norm + # canonical + if compare_canonical is True: + dbi_eval = deepcopy(dbi_object) + dbi_eval.mode = DoubleBracketGeneratorType.canonical + if step is None: + step_best = dbi_eval.choose_step(scheduling=scheduling, **kwargs) + else: + step_best = step + dbi_eval(step=step_best) + optimal_steps[-1] = step_best + norms_off_diagonal_restriction[-1] = dbi_eval.off_diagonal_norm + # find best d + idx_max_loss = np.argmin(norms_off_diagonal_restriction) + flip = flip_list[idx_max_loss] + step_optimal = optimal_steps[idx_max_loss] + dbi_eval = deepcopy(dbi_object) + if idx_max_loss == len(d_list) and compare_canonical is True: + # canonical + dbi_eval(step=step_optimal, mode=DoubleBracketGeneratorType.canonical) + + else: + d_optimal = flip * d_list[idx_max_loss] + dbi_eval(step=step_optimal, d=d_optimal) + return dbi_eval, idx_max_loss, step_optimal, flip + + +def gradient_numerical( + dbi_object: DoubleBracketIteration, + d_params: list, + parameterization: ParameterizationTypes, + s: float = 1e-2, + delta: float = 1e-3, + backend=None, + **kwargs, +): + r""" + Gradient of the DBI with respect to the parametrization of D. A simple finite difference is used to calculate the gradient. + + Args: + dbi_object (DoubleBracketIteration): DoubleBracketIteration object. + d_params (np.array): Parameters for the ansatz (note that the dimension must be 2**nqubits for full ansazt and nqubits for Pauli ansatz). + s (float): A short flow duration for finding the numerical gradient. + delta (float): Step size for numerical gradient. + Returns: + grad (np.array): Gradient of the D operator. + """ + backend = _check_backend(backend) + nqubits = dbi_object.nqubits + grad = np.zeros(len(d_params)) + d = params_to_diagonal_operator( + d_params, nqubits, parameterization=parameterization, **kwargs + ) + for i in range(len(d_params)): + params_new = backend.cast(d_params, copy=True) + params_new[i] += delta + d_new = params_to_diagonal_operator( + params_new, nqubits, parameterization=parameterization, **kwargs + ) + # find the increment of a very small step + grad[i] = (dbi_object.loss(s, d_new) - dbi_object.loss(s, d)) / delta + return grad + + +def gradient_descent( + dbi_object: DoubleBracketIteration, + iterations: int, + d_params: list, + parameterization: ParameterizationTypes, + pauli_operator_dict: dict = None, + pauli_parameterization_order: int = 1, + normalize: bool = False, + lr_min: float = 1e-5, + lr_max: float = 1, + max_evals: int = 100, + space: callable = None, + optimizer: callable = hyperopt.tpe, + verbose: bool = False, + backend=None, +): + r"""Numerical gradient descent method for variating diagonal operator in each double bracket rotation. + + Args: + dbi_object (DoubleBracketIteration): the target double bracket object. + iterations (int): number of double bracket rotations. + d_params (list): the parameters for the initial diagonal operator. + parameterization (ParameterizationTypes): the parameterization method for diagonal operator. + Options include pauli and computational. + pauli_operator_dict (dictionary, optional): dictionary of "name": Pauli-operator for Pauli-based parameterization type. + Defaults to None. + pauli_parameterization_order (int, optional): the order of parameterization or locality in Pauli basis. Defaults to 1. + normalize (bool, optional): option to normalize the diagonal operator. Defaults to False. + lr_min (float, optional): the minimal gradient step. Defaults to 1e-5. + lr_max (float, optional): the maximal gradient step. Defaults to 1. + max_evals (int, optional): maximum number of evaluations for `lr` using `hyperopt`. Defaults to 100. + space (callable, optional): evalutation space for `hyperopt`. Defaults to None. + optimizer (callable, optional): optimizer option for `hyperopt`. Defaults to `hyperopt.tpe`. + verbose (bool, optional): option for printing `hyperopt` process. Defaults to False. + + Returns: + loss_hist (list): list of history losses of `dbi_object` throughout the double bracket rotations. + d_params_hist (list): list of history of `d` parameters after gradient descent. + s_hist (list): list of history of optimal `s` found. + Example: + from qibo import set_backend + from qibo.hamiltonians import Hamiltonian + from qibo.models.dbi.double_bracket import * + from qibo.models.dbi.utils import * + from qibo.models.dbi.utils_dbr_strategies import gradient_descent + from qibo.quantum_info import random_hermitian + + nqubits = 3 + NSTEPS = 5 + set_backend("numpy") + h0 = random_hermitian(2**nqubits) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0), + mode=DoubleBracketGeneratorType.single_commutator, + scheduling=DoubleBracketScheduling.hyperopt, + cost=DoubleBracketCostFunction.off_diagonal_norm, + ) + initial_off_diagonal_norm = dbi.off_diagonal_norm + pauli_operator_dict = generate_pauli_operator_dict( + nqubits, parameterization_order=1 + ) + pauli_operators = list(pauli_operator_dict.values()) + # let initial d be approximation of $\Delta(H) + d_coef_pauli = decompose_into_Pauli_basis( + dbi.diagonal_h_matrix, pauli_operators=pauli_operators + ) + d_pauli = sum([d_coef_pauli[i] * pauli_operators[i] for i in range(nqubits)]) + loss_hist_pauli, d_params_hist_pauli, s_hist_pauli = gradient_descent( + dbi, + NSTEPS, + d_coef_pauli, + ParameterizationTypes.pauli, + pauli_operator_dict=pauli_operator_dict, + ) + """ + backend = _check_backend(backend) + + nqubits = dbi_object.nqubits + # TODO: write tests where this condition applies + if ( + parameterization is ParameterizationTypes.pauli and pauli_operator_dict is None + ): # pragma: no cover + pauli_operator_dict = generate_pauli_operator_dict( + nqubits=nqubits, parameterization_order=pauli_parameterization_order + ) + d = params_to_diagonal_operator( + d_params, + nqubits, + parameterization=parameterization, + pauli_operator_dict=pauli_operator_dict, + normalize=normalize, + ) + loss_hist = [dbi_object.loss(0.0, d=d)] + d_params_hist = [d_params] + s_hist = [0] + # TODO: write tests where this condition applies + if ( + parameterization is ParameterizationTypes.pauli and pauli_operator_dict is None + ): # pragma: no cover + pauli_operator_dict = generate_pauli_operator_dict( + nqubits=nqubits, + parameterization_order=pauli_parameterization_order, + backend=backend, + ) + # first step + s = dbi_object.choose_step(d=d) + dbi_object(step=s, d=d) + for _ in range(iterations): + grad = gradient_numerical( + dbi_object, + d_params, + parameterization, + pauli_operator_dict=pauli_operator_dict, + pauli_parameterization_order=pauli_parameterization_order, + normalize=normalize, + backend=backend, + ) + + # set up hyperopt to find optimal lr + def func_loss_to_lr(lr): + d_params_eval = [d_params[j] - grad[j] * lr for j in range(len(grad))] + d_eval = params_to_diagonal_operator( + d_params_eval, + nqubits, + parameterization=parameterization, + pauli_operator_dict=pauli_operator_dict, + normalize=normalize, + ) + return dbi_object.loss(step=s, d=d_eval) + + if space is None: + space = hyperopt.hp.loguniform("lr", np.log(lr_min), np.log(lr_max)) + + best = hyperopt.fmin( + fn=func_loss_to_lr, + space=space, + algo=optimizer.suggest, + max_evals=max_evals, + verbose=verbose, + ) + lr = best["lr"] + + d_params = [d_params[j] - grad[j] * lr for j in range(len(grad))] + d = params_to_diagonal_operator( + d_params, + nqubits, + parameterization=parameterization, + pauli_operator_dict=pauli_operator_dict, + normalize=normalize, + ) + s = dbi_object.choose_step(d=d) + dbi_object(step=s, d=d) + + # record history + loss_hist.append(dbi_object.loss(0.0, d=d)) + d_params_hist.append(d_params) + s_hist.append(s) + return loss_hist, d_params_hist, s_hist diff --git a/src/qibo/models/dbi/utils_scheduling.py b/src/qibo/models/dbi/utils_scheduling.py new file mode 100644 index 0000000000..130cd88f30 --- /dev/null +++ b/src/qibo/models/dbi/utils_scheduling.py @@ -0,0 +1,212 @@ +import math +from functools import partial +from typing import Optional + +import hyperopt +import numpy as np + +error = 1e-3 + + +def grid_search_step( + dbi_object, + step_min: float = 1e-5, + step_max: float = 1, + num_evals: int = 100, + space: Optional[np.array] = None, + d: Optional[np.array] = None, +): + """ + Greedy optimization of the iteration step. + + Args: + step_min: lower bound of the search grid; + step_max: upper bound of the search grid; + mnum_evals: number of iterations between step_min and step_max; + d: diagonal operator for generating double-bracket iterations. + + Returns: + (float): optimized best iteration step (minimizing off-diagonal norm). + """ + if space is None: + space = np.linspace(step_min, step_max, num_evals) + + if d is None: + d = dbi_object.diagonal_h_matrix + + loss_list = [dbi_object.loss(step, d=d) for step in space] + + idx_max_loss = np.argmin(loss_list) + return space[idx_max_loss] + + +def hyperopt_step( + dbi_object, + step_min: float = 1e-5, + step_max: float = 1, + max_evals: int = 100, + space: callable = None, + optimizer: callable = None, + look_ahead: int = 1, + d: Optional[np.array] = None, +): + """ + Optimize iteration step using hyperopt. + + Args: + step_min: lower bound of the search grid; + step_max: upper bound of the search grid; + max_evals: maximum number of iterations done by the hyperoptimizer; + space: see hyperopt.hp possibilities; + optimizer: see hyperopt algorithms; + look_ahead: number of iteration steps to compute the loss function; + d: diagonal operator for generating double-bracket iterations. + + Returns: + (float): optimized best iteration step (minimizing loss function). + """ + if space is None: + space = hyperopt.hp.uniform + if optimizer is None: + optimizer = hyperopt.tpe + if d is None: + d = dbi_object.diagonal_h_matrix + + space = space("step", step_min, step_max) + + best = hyperopt.fmin( + fn=partial(dbi_object.loss, d=d, look_ahead=look_ahead), + space=space, + algo=optimizer.suggest, + max_evals=max_evals, + show_progressbar=False, + ) + return best["step"] + + +def polynomial_step( + dbi_object, + n: int = 2, + n_max: int = 5, + d: np.array = None, + coef: Optional[list] = None, + cost: Optional[str] = None, +): + r""" + Optimizes iteration step by solving the n_th order polynomial expansion of the loss function. + e.g. $n=2$: $2\Trace(\sigma(\Gamma_1 + s\Gamma_2 + s^2/2\Gamma_3)\sigma(\Gamma_0 + s\Gamma_1 + s^2/2\Gamma_2)) + Args: + n (int, optional): the order to which the loss function is expanded. Defaults to 4. + n_max (int, optional): maximum order allowed for recurring calls of `polynomial_step`. Defaults to 5. + d (np.array, optional): diagonal operator, default as $\delta(H)$. + backup_scheduling (`DoubleBracketScheduling`): the scheduling method to use in case no real positive roots are found. + """ + if cost is None: + cost = dbi_object.cost + + if d is None: + d = dbi_object.diagonal_h_matrix + + if n > n_max: + raise ValueError( + "No solution can be found with polynomial approximation. Increase `n_max` or use other scheduling methods." + ) + if coef is None: + coef = dbi_object.cost_expansion(d=d, n=n) + roots = np.roots(coef) + real_positive_roots = [ + np.real(root) for root in roots if np.imag(root) < 1e-3 and np.real(root) > 0 + ] + # solution exists, return minimum s + if len(real_positive_roots) > 0: + losses = [dbi_object.loss(step=root, d=d) for root in real_positive_roots] + return real_positive_roots[losses.index(min(losses))] + # solution does not exist, return None + else: + return None + + +def simulated_annealing_step( + dbi_object, + d: Optional[np.array] = None, + initial_s=None, + step_min=1e-5, + step_max=1, + s_jump_range=None, + s_jump_range_divident=5, + initial_temp=1, + cooling_rate=0.85, + min_temp=1e-5, + max_iter=200, +): + """ + Perform a single step of simulated annealing optimization. + + Parameters: + dbi_object: DBI object + The object representing the problem to be optimized. + d: Optional[np.array], optional + The diagonal matrix 'd' used in optimization. If None, it uses the diagonal + matrix 'diagonal_h_matrix' from dbi_object. + initial_s: float or None, optional + Initial value for 's', the step size. If None, it is initialized using + polynomial_step function with 'n=4'. If 'polynomial_step' returns None, + 'initial_s' is set to 'step_min'. + step_min: float, optional + Minimum value for the step size 's'. + step_max: float, optional + Maximum value for the step size 's'. + s_jump_range: float or None, optional + Range for the random jump in step size. If None, it's calculated based on + 'step_min', 'step_max', and 's_jump_range_divident'. + s_jump_range_divident: int, optional + Dividend to determine the range for random jump in step size. + initial_temp: float, optional + Initial temperature for simulated annealing. + cooling_rate: float, optional + Rate at which temperature decreases in simulated annealing. + min_temp: float, optional + Minimum temperature threshold for termination of simulated annealing. + max_iter: int, optional + Maximum number of iterations for simulated annealing. + + Returns: + float: + The optimized step size 's'. + """ + + if d is None: + d = dbi_object.diagonal_h_matrix + if initial_s is None: + initial_s = polynomial_step(dbi_object=dbi_object, d=d, n=4) + # TODO: implement test to catch this if statement + if initial_s is None: # pragma: no cover + initial_s = step_min + if s_jump_range is None: + s_jump_range = (step_max - step_min) / s_jump_range_divident + current_s = initial_s + current_loss = dbi_object.loss(d=d, step=current_s) + temp = initial_temp + + for _ in range(max_iter): + candidate_s = max( + step_min, + min( + current_s + np.random.uniform(-1 * s_jump_range, s_jump_range), step_max + ), + ) + candidate_loss = dbi_object.loss(d=d, step=candidate_s) + + # Calculate change in loss + delta_loss = candidate_loss - current_loss + + # Determine if the candidate solution is an improvement + if delta_loss < 0 or np.random.rand() < math.exp(-delta_loss / temp): + current_s = candidate_s + current_loss = candidate_loss + # Cool down + temp *= cooling_rate + if temp < min_temp or current_s > step_max or current_s < step_min: + break + + return current_s diff --git a/tests/test_models_dbi.py b/tests/test_models_dbi.py index 4a2133f33c..3ce40e6559 100644 --- a/tests/test_models_dbi.py +++ b/tests/test_models_dbi.py @@ -3,21 +3,31 @@ import numpy as np import pytest +from qibo import hamiltonians, set_backend from qibo.hamiltonians import Hamiltonian from qibo.models.dbi.double_bracket import ( + DoubleBracketCostFunction, DoubleBracketGeneratorType, DoubleBracketIteration, + DoubleBracketScheduling, ) +from qibo.models.dbi.utils import * +from qibo.models.dbi.utils_dbr_strategies import ( + gradient_descent, + select_best_dbr_generator, +) +from qibo.models.dbi.utils_scheduling import polynomial_step from qibo.quantum_info import random_hermitian -NSTEPS = 50 -SEED = 10 +NSTEPS = 3 +seed = 10 """Number of steps for evolution.""" @pytest.mark.parametrize("nqubits", [1, 2]) def test_double_bracket_iteration_canonical(backend, nqubits): - h0 = random_hermitian(2**nqubits, backend=backend) + """Check default (canonical) mode.""" + h0 = random_hermitian(2**nqubits, backend=backend, seed=seed) dbi = DoubleBracketIteration( Hamiltonian(nqubits, h0, backend=backend), mode=DoubleBracketGeneratorType.canonical, @@ -31,7 +41,8 @@ def test_double_bracket_iteration_canonical(backend, nqubits): @pytest.mark.parametrize("nqubits", [1, 2]) def test_double_bracket_iteration_group_commutator(backend, nqubits): - h0 = random_hermitian(2**nqubits, backend=backend, seed=SEED) + """Check group commutator mode.""" + h0 = random_hermitian(2**nqubits, backend=backend, seed=seed) d = backend.cast(np.diag(np.diag(backend.to_numpy(h0)))) dbi = DoubleBracketIteration( Hamiltonian(nqubits, h0, backend=backend), @@ -41,7 +52,6 @@ def test_double_bracket_iteration_group_commutator(backend, nqubits): # test first iteration with default d dbi(mode=DoubleBracketGeneratorType.group_commutator, step=0.01) - for _ in range(NSTEPS): dbi(step=0.01, d=d) @@ -49,8 +59,8 @@ def test_double_bracket_iteration_group_commutator(backend, nqubits): @pytest.mark.parametrize("nqubits", [1, 2]) -def test_double_bracket_iteration_group_commutator_3(backend, nqubits): - """Check group commutator mode.""" +def test_double_bracket_iteration_group_commutator_3rd_order(backend, nqubits): + """Check 3rd order group commutator mode.""" h0 = random_hermitian(2**nqubits, backend=backend, seed=seed) d = backend.cast(np.diag(np.diag(backend.to_numpy(h0)))) dbi = DoubleBracketIteration( @@ -69,7 +79,8 @@ def test_double_bracket_iteration_group_commutator_3(backend, nqubits): @pytest.mark.parametrize("nqubits", [1, 2]) def test_double_bracket_iteration_single_commutator(backend, nqubits): - h0 = random_hermitian(2**nqubits, backend=backend) + """Check single commutator mode.""" + h0 = random_hermitian(2**nqubits, backend=backend, seed=seed) d = backend.cast(np.diag(np.diag(backend.to_numpy(h0)))) dbi = DoubleBracketIteration( Hamiltonian(nqubits, h0, backend=backend), @@ -86,49 +97,190 @@ def test_double_bracket_iteration_single_commutator(backend, nqubits): assert initial_off_diagonal_norm > dbi.off_diagonal_norm -@pytest.mark.parametrize("nqubits", [3, 4]) -def test_hyperopt_step(backend, nqubits): - h0 = random_hermitian(2**nqubits, backend=backend) - d = backend.cast(np.diag(np.diag(backend.to_numpy(h0)))) - dbi = DoubleBracketIteration(Hamiltonian(nqubits, h0, backend=backend)) +@pytest.mark.parametrize("nqubits", [2, 3]) +@pytest.mark.parametrize( + "scheduling", + [ + DoubleBracketScheduling.grid_search, + DoubleBracketScheduling.hyperopt, + DoubleBracketScheduling.simulated_annealing, + ], +) +def test_variational_scheduling(backend, nqubits, scheduling): + """Check schduling options.""" + h = 2 + # define the hamiltonian + h0 = hamiltonians.TFIM(nqubits=nqubits, h=h) + dbi = DoubleBracketIteration(h0, scheduling=scheduling) # find initial best step with look_ahead = 1 - initial_step = 0.01 - delta = 0.02 + initial_off_diagonal_norm = dbi.off_diagonal_norm + for _ in range(NSTEPS): + step = dbi.choose_step() + dbi(step=step) + assert initial_off_diagonal_norm > dbi.off_diagonal_norm + - step = dbi.hyperopt_step( - step_min=initial_step - delta, step_max=initial_step + delta, max_evals=10 +@pytest.mark.parametrize( + "cost", + [ + DoubleBracketCostFunction.off_diagonal_norm, + DoubleBracketCostFunction.least_squares, + ], +) +def test_polynomial_cost_function(backend, cost): + nqubits = 2 + h0 = random_hermitian(2**nqubits, backend=backend, seed=seed) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0, backend=backend), + mode=DoubleBracketGeneratorType.single_commutator, + cost=cost, + scheduling=DoubleBracketScheduling.polynomial_approximation, ) + initial_off_diagonal_norm = dbi.off_diagonal_norm + for i in range(NSTEPS): + s = dbi.choose_step(d=dbi.diagonal_h_matrix, n=5) + dbi(step=s, d=dbi.off_diag_h) + assert initial_off_diagonal_norm > dbi.off_diagonal_norm - assert step != initial_step - # evolve following the optimized first step - for generator in DoubleBracketGeneratorType: - dbi(mode=generator, step=step, d=d) +def test_polynomial_energy_fluctuation(backend): + nqubits = 4 + h0 = random_hermitian(2**nqubits, seed=seed, backend=backend) + state = np.zeros(2**nqubits) + state[0] = 1 + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0, backend=backend), + mode=DoubleBracketGeneratorType.single_commutator, + cost=DoubleBracketCostFunction.energy_fluctuation, + scheduling=DoubleBracketScheduling.polynomial_approximation, + ref_state=state, + ) + for i in range(NSTEPS): + s = dbi.choose_step(d=dbi.diagonal_h_matrix, n=5) + dbi(step=s, d=dbi.diagonal_h_matrix) + assert dbi.energy_fluctuation(state=state) < dbi.h0.energy_fluctuation(state=state) - # find the following step size with look_ahead - look_ahead = 3 - step = dbi.hyperopt_step( - step_min=initial_step - delta, - step_max=initial_step + delta, - max_evals=10, - look_ahead=look_ahead, +@pytest.mark.parametrize("nqubits", [5, 6]) +def test_polynomial_fail_cases(backend, nqubits): + h0 = random_hermitian(2**nqubits, backend=backend, seed=seed) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0, backend=backend), + mode=DoubleBracketGeneratorType.single_commutator, + scheduling=DoubleBracketScheduling.polynomial_approximation, ) + with pytest.raises(ValueError): + polynomial_step(dbi, n=2, n_max=1) + assert polynomial_step(dbi, n=1) == None - # evolve following the optimized first step - for gentype in range(look_ahead): - dbi(mode=DoubleBracketGeneratorType(gentype + 1), step=step, d=d) +def test_least_squares(backend): + """Check least squares cost function.""" + nqubits = 4 + h0 = random_hermitian(2**nqubits, backend=backend, seed=seed) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0, backend=backend), + cost=DoubleBracketCostFunction.least_squares, + ) + d = np.diag(np.linspace(1, 2**nqubits, 2**nqubits)) / 2**nqubits + initial_potential = dbi.least_squares(d=d) + step = dbi.choose_step(d=d) + dbi(d=d, step=step) + assert dbi.least_squares(d=d) < initial_potential -def test_energy_fluctuations(backend): - h0 = np.array([[1, 0], [0, -1]]) - h0 = backend.cast(h0, dtype=backend.dtype) - state = np.array([1, 0]) - state = backend.cast(state, dtype=backend.dtype) +@pytest.mark.parametrize("compare_canonical", [True, False]) +@pytest.mark.parametrize("step", [None, 1e-3]) +@pytest.mark.parametrize("nqubits", [2, 3]) +def test_select_best_dbr_generator(backend, nqubits, step, compare_canonical): + h0 = random_hermitian(2**nqubits, backend=backend, seed=seed) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0, backend=backend), + mode=DoubleBracketGeneratorType.single_commutator, + ) + initial_off_diagonal_norm = dbi.off_diagonal_norm + generate_local_Z = generate_Z_operators(nqubits, backend=backend) + Z_ops = list(generate_local_Z.values()) + for _ in range(NSTEPS): + dbi, idx, step, flip_sign = select_best_dbr_generator( + dbi, + Z_ops, + compare_canonical=compare_canonical, + step=step, + ) + assert dbi.off_diagonal_norm < initial_off_diagonal_norm + + +@pytest.mark.parametrize("step", [None, 1e-3]) +def test_params_to_diagonal_operator(backend, step): + nqubits = 2 + pauli_operator_dict = generate_pauli_operator_dict( + nqubits, parameterization_order=1, backend=backend + ) + params = [1, 2, 3] + operator_pauli = sum( + [params[i] * list(pauli_operator_dict.values())[i] for i in range(nqubits)] + ) + backend.assert_allclose( + operator_pauli, + params_to_diagonal_operator( + params, + nqubits=nqubits, + parameterization=ParameterizationTypes.pauli, + pauli_operator_dict=pauli_operator_dict, + ), + ) + operator_element = params_to_diagonal_operator( + params, + nqubits=nqubits, + parameterization=ParameterizationTypes.computational, + ) + for i in range(len(params)): + backend.assert_allclose( + backend.cast(backend.to_numpy(operator_element).diagonal())[i], params[i] + ) - dbi = DoubleBracketIteration(Hamiltonian(1, matrix=h0, backend=backend)) - energy_fluctuation = dbi.energy_fluctuation(state=state) - assert energy_fluctuation == 0 +@pytest.mark.parametrize("order", [1, 2]) +def test_gradient_descent(backend, order): + nqubits = 2 + h0 = random_hermitian(2**nqubits, seed=seed, backend=backend) + dbi = DoubleBracketIteration( + Hamiltonian(nqubits, h0, backend=backend), + mode=DoubleBracketGeneratorType.single_commutator, + scheduling=DoubleBracketScheduling.hyperopt, + cost=DoubleBracketCostFunction.off_diagonal_norm, + ) + initial_off_diagonal_norm = dbi.off_diagonal_norm + pauli_operator_dict = generate_pauli_operator_dict( + nqubits, + parameterization_order=order, + backend=backend, + ) + pauli_operators = list(pauli_operator_dict.values()) + # let initial d be approximation of $\Delta(H) + d_coef_pauli = decompose_into_pauli_basis( + dbi.diagonal_h_matrix, pauli_operators=pauli_operators + ) + d_pauli = sum([d_coef_pauli[i] * pauli_operators[i] for i in range(nqubits)]) + loss_hist_pauli, d_params_hist_pauli, s_hist_pauli = gradient_descent( + dbi, + NSTEPS, + d_coef_pauli, + ParameterizationTypes.pauli, + pauli_operator_dict=pauli_operator_dict, + pauli_parameterization_order=order, + ) + assert loss_hist_pauli[-1] < initial_off_diagonal_norm + + # computational basis + d_coef_computational_partial = backend.cast(backend.to_numpy(d_pauli).diagonal()) + ( + loss_hist_computational_partial, + _, + _, + ) = gradient_descent( + dbi, NSTEPS, d_coef_computational_partial, ParameterizationTypes.computational + ) + assert loss_hist_computational_partial[-1] < initial_off_diagonal_norm diff --git a/tests/test_models_dbi_utils.py b/tests/test_models_dbi_utils.py deleted file mode 100644 index 9f55171547..0000000000 --- a/tests/test_models_dbi_utils.py +++ /dev/null @@ -1,51 +0,0 @@ -""""Testing utils for DoubleBracketIteration model""" - -import numpy as np -import pytest - -from qibo.hamiltonians import Hamiltonian -from qibo.models.dbi.double_bracket import ( - DoubleBracketGeneratorType, - DoubleBracketIteration, -) -from qibo.models.dbi.utils import * -from qibo.quantum_info import random_hermitian - -NSTEPS = 5 -"""Number of steps for evolution.""" - - -@pytest.mark.parametrize("nqubits", [1, 2]) -def test_generate_Z_operators(backend, nqubits): - h0 = random_hermitian(2**nqubits, backend=backend) - dbi = DoubleBracketIteration( - Hamiltonian(nqubits=nqubits, matrix=h0, backend=backend) - ) - generate_Z = generate_Z_operators(nqubits, backend=backend) - Z_ops = list(generate_Z.values()) - - delta_h0 = dbi.diagonal_h_matrix - dephasing_channel = (sum([Z_op @ h0 @ Z_op for Z_op in Z_ops]) + h0) / 2**nqubits - norm_diff = np.linalg.norm(delta_h0 - dephasing_channel) - - assert norm_diff < 1e-3 - - -@pytest.mark.parametrize("nqubits", [1, 2]) -@pytest.mark.parametrize("step", [0.1, 0.2]) -def test_select_best_dbr_generator(backend, nqubits, step): - h0 = random_hermitian(2**nqubits, seed=1, backend=backend) - dbi = DoubleBracketIteration( - Hamiltonian(nqubits, h0, backend=backend), - mode=DoubleBracketGeneratorType.single_commutator, - ) - generate_Z = generate_Z_operators(nqubits, backend=backend) - Z_ops = list(generate_Z.values()) - initial_off_diagonal_norm = dbi.off_diagonal_norm - - for _ in range(NSTEPS): - dbi, idx, step_optimize, flip = select_best_dbr_generator( - dbi, Z_ops, step=step, compare_canonical=True, max_evals=5 - ) - - assert initial_off_diagonal_norm > dbi.off_diagonal_norm