diff --git a/BIBLIOGRAPHY.md b/BIBLIOGRAPHY.md index d683f291..ef2a7b8a 100644 --- a/BIBLIOGRAPHY.md +++ b/BIBLIOGRAPHY.md @@ -16,7 +16,8 @@ The code in this project is based on ideas from the following publications: - **[FowlerWright2022]** Fowler-Wright et al., *Efficient Many-Body Non-Markovian Dynamics of Organic Polaritons*, [Phys. Rev. Lett. 129, 173001](https://doi.org/10.1103/PhysRevLett.129.173001) (2022). - **[Fux2023]** Fux et al., *Tensor network simulation of chains of non-Markovian open quantum systems*, [Phys. Rev. Research 5, 033078 ](https://doi.org/10.1103/PhysRevResearch.5.033078}) (2023). - **[Butler2024]** Butler et al., *Optimizing Performance of Quantum Operations with Non-Markovian Decoherence: The Tortoise or the Hare?*, [Phys. Rev. Lett. 132, 060401 ](https://doi.org/10.1103/PhysRevLett.132.060401}) (2024). - +- **[Link2024]** Link et al., *Open Quantum System Dynamics from Infinite Tensor Network Contraction*, `Phys. Rev. Lett. 132, + 200403 `__ (2024). BibTeX: ------- @@ -142,6 +143,21 @@ BibTeX: url = {https://doi.org/10.1103/PhysRevLett.123.240602} } +@article{Link2024, + title = {Open Quantum System Dynamics from Infinite Tensor Network Contraction}, + author = {Link, Valentin and Tu, Hong-Hao and Strunz, Walter T.}, + journal = {Phys. Rev. Lett.}, + volume = {132}, + issue = {20}, + pages = {200403}, + numpages = {7}, + year = {2024}, + month = {May}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevLett.132.200403}, + url = {https://link.aps.org/doi/10.1103/PhysRevLett.132.200403} +} + @misc{OQuPy-code-v0.5.0, author={{The TEMPO collaboration}}, title={{OQuPy: A Python package to efficiently simulate non-Markovian open quantum systems with process tensors.}}, diff --git a/HOW_TO_CITE.md b/HOW_TO_CITE.md index f082d91b..c34bf3fa 100644 --- a/HOW_TO_CITE.md +++ b/HOW_TO_CITE.md @@ -20,7 +20,7 @@ Also, please consider citing: - Chains (PT-TEBD): **[Fux2023]** - Mean-Field TEMPO: **[FowlerWright2022]** - Gibbs TEMPO: **[Chiu2022]** - +- Time translational invariant TEMPO: **[Link2024]** BibTeX ------ diff --git a/README.md b/README.md index 85e106d2..17db6c00 100644 --- a/README.md +++ b/README.md @@ -41,7 +41,7 @@ Furthermore, OQuPy implements methods to ... - **[9]** Butler et al., [Phys. Rev. Lett. 132, 060401 ](https://doi.org/10.1103/PhysRevLett.132.060401) (2024). - **[10]** Gribben et al., [Quantum, 6, 847](https://doi.org/10.22331/q-2022-10-25-847) (2022). - **[11]** Chiu et al., [Phys. Rev. A 106, 012204](https://doi.org/10.1103/PhysRevA.106.012204) (2022). - +- **[12]** Link et al., [Phys. Rev. Lett. 132, 200403](https://doi.org/10.1103/PhysRevLett.132.200403) (2024). ------------------------------------------------------------------------------- diff --git a/docs/index.rst b/docs/index.rst index 16e2eaf2..029c19e6 100644 --- a/docs/index.rst +++ b/docs/index.rst @@ -61,7 +61,8 @@ Furthermore, OQuPy implements methods to ... 847 `__ (2022). - **[11]** Chiu et al., `Phys. Rev. A 106, 012204 `__ (2022). - +- **[12]** Link et al., `Phys. Rev. Lett. 132, + 200403 `__ (2024). .. |binder-tutorial| image:: https://mybinder.org/badge_logo.svg :target: https://mybinder.org/v2/gh/tempoCollaboration/OQuPy/main?filepath=tutorials%2Fquickstart.ipynb diff --git a/docs/pages/bibliography.rst b/docs/pages/bibliography.rst index d4c30356..fdc28a36 100644 --- a/docs/pages/bibliography.rst +++ b/docs/pages/bibliography.rst @@ -45,6 +45,9 @@ publications: Operations with Non-Markovian Decoherence: The Tortoise or the Hare?*, `Phys. Rev. Lett. 132, 060401 `__ (2024). +- **[Link2024]** Link et al., *Open Quantum System Dynamics from Infinite + Tensor Network Contraction*, `Phys. Rev. Lett. 132, + 200403 `__ (2024). .. _bibtex: @@ -173,6 +176,21 @@ BibTeX: url = {https://doi.org/10.1103/PhysRevLett.123.240602} } + @article{Link2024, + title = {Open Quantum System Dynamics from Infinite Tensor Network Contraction}, + author = {Link, Valentin and Tu, Hong-Hao and Strunz, Walter T.}, + journal = {Phys. Rev. Lett.}, + volume = {132}, + issue = {20}, + pages = {200403}, + numpages = {7}, + year = {2024}, + month = {May}, + publisher = {American Physical Society}, + doi = {10.1103/PhysRevLett.132.200403}, + url = {https://link.aps.org/doi/10.1103/PhysRevLett.132.200403} + } + @misc{OQuPy-code-v0.5.0, author={{The TEMPO collaboration}}, title={{OQuPy: A Python package to efficiently simulate non-Markovian open quantum systems with process tensors.}}, diff --git a/docs/pages/how_to_cite.rst b/docs/pages/how_to_cite.rst index a00f1063..fd68401c 100644 --- a/docs/pages/how_to_cite.rst +++ b/docs/pages/how_to_cite.rst @@ -20,6 +20,7 @@ Also, please consider citing: - Chains (PT-TEBD): **[Fux2023]** - Mean-Field TEMPO: **[FowlerWright2022]** - Gibbs TEMPO: **[Chiu2022]** +- Time translational invariant TEMPO: **[Link2024]** BibTeX ------ diff --git a/examples/fidelity-gradient-tti.ipynb b/examples/fidelity-gradient-tti.ipynb new file mode 100644 index 00000000..a8d149a4 --- /dev/null +++ b/examples/fidelity-gradient-tti.ipynb @@ -0,0 +1,291 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# %load examples/fidelity-gradient-tti.py\n", + "#!/usr/bin/env python\n", + "\n", + "# Computation takes roughly 10 minutes.\n", + "\n", + "import sys\n", + "sys.path.insert(0,'..')\n", + "\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "import oqupy\n", + "import oqupy.operators as op\n", + "\n", + "from oqupy.tti_tempo import TTITempo\n", + "from scipy.optimize import minimize\n", + "\n", + "# --- Parameters --------------------------------------------------------------\n", + "\n", + "# parameters used in https://arxiv.org/abs/2303.16002\n", + "\n", + "# -- time steps --\n", + "end_time = 2.5\n", + "num_steps = 50\n", + "dt = end_time / num_steps\n", + "\n", + "num_steps = 200\n", + "\n", + "# -- bath --\n", + "alpha = 0.126 \n", + "omega_cutoff = 3.04 \n", + "temperature = 5 * 0.1309 \n", + "pt_dkmax = 60 \n", + "pt_epsrel = 10**(-7) \n", + "\n", + "# -- initial and target state --\n", + "initial_state = op.spin_dm('x-')\n", + "target_derivative = op.spin_dm('x+').T\n", + "\n", + "# -- initial parameter guess, defined over half-timesteps --\n", + "y0 = np.zeros(2*num_steps)\n", + "z0 = np.ones(2*num_steps) * (np.pi) / (num_steps*dt)\n", + "x0 = np.zeros(2*num_steps)\n", + "\n", + "parameters = np.vstack((x0,y0,z0)).T\n", + "num_params = parameters.shape[1]\n", + "print(parameters.shape)\n", + "\n", + "# --- Compute process tensors -------------------------------------------------\n", + "\n", + "correlations = oqupy.PowerLawSD(\n", + " alpha=alpha,\n", + " zeta=3,\n", + " cutoff=omega_cutoff,\n", + " cutoff_type='gaussian',\n", + " temperature=temperature)\n", + "bath = oqupy.Bath(0.5* op.sigma('z'), correlations)\n", + "pt_tempo_parameters = oqupy.TempoParameters(\n", + " dt=dt,\n", + " epsrel=pt_epsrel,\n", + " dkmax=pt_dkmax)\n", + "\n", + "pt=TTITempo(bath,start_time=0.0,parameters=pt_tempo_parameters)\n", + "process_tensor=pt.get_process_tensor()\n", + "\n", + "#process_tensor = oqupy.pt_tempo_compute(\n", + "# bath=bath,\n", + "# start_time=0.0,\n", + "# end_time=pt_end_time,\n", + "# parameters=pt_tempo_parameters)\n", + "\n", + "# --- Define parameterized system ----------------------------------------------\n", + "\n", + "hs_dim = 2\n", + "\n", + "def hamiltonian(x, y, z):\n", + " h = np.zeros((2,2),dtype='complex128')\n", + " for var, var_name in zip([x,y,z], ['x', 'y', 'z']):\n", + " h += 0.5* var * op.sigma(var_name)\n", + " return h\n", + "\n", + "parametrized_system = oqupy.ParameterizedSystem(hamiltonian)\n", + "\n", + "# --- Compute fidelity, dynamics, and fidelity gradient -----------------------\n", + "\n", + "from oqupy.gradient import state_gradient\n", + "\n", + "fidelity_dict = state_gradient(\n", + " system=parametrized_system,\n", + " initial_state=initial_state,\n", + " target_derivative=target_derivative,\n", + " process_tensors=[process_tensor],\n", + " parameters=parameters,\n", + " num_steps=num_steps)\n", + "\n", + "final_state = fidelity_dict['dynamics'].states[-1]\n", + "v_final_state = np.reshape(final_state,hs_dim**2)\n", + "v_target_derivative = np.reshape(target_derivative.T,hs_dim**2)\n", + "fidelity = v_target_derivative @ v_final_state\n", + "\n", + "print(f\"the fidelity is {fidelity.real}\")\n", + "#print(f\"the fidelity gradient is {fidelity_dict['gradient'].real}\")\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "t, s_x = fidelity_dict['dynamics'].expectations(op.sigma(\"x\"), real=True)\n", + "t, s_y = fidelity_dict['dynamics'].expectations(op.sigma(\"y\"), real=True)\n", + "t, s_z = fidelity_dict['dynamics'].expectations(op.sigma(\"z\"), real=True)\n", + "bloch_length = np.sqrt(s_x**2 +s_y**2 + s_z**2)\n", + "\n", + "plt.title(\"Pre-optimisation evolution of components\")\n", + "plt.plot(t,s_x,label='x')\n", + "plt.plot(t,s_y,label='y')\n", + "plt.plot(t,s_z,label='z')\n", + "plt.plot(t,bloch_length,label=r'$|\\mathbf{\\sigma}|$')\n", + "\n", + "plt.ylabel(r\"$\\langle \\sigma \\rangle$\",rotation=0,fontsize=16)\n", + "plt.xlabel(\"t\")\n", + "\n", + "plt.legend()\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "# ----------- Optimisation of control parameters w.r.t. infidelity ---------------\n", + "\n", + "def infidandgrad(paras):\n", + " \"\"\"\"\"\n", + " Take a numpy array [hx0, hz0, hx1, hz1, ...] over full timesteps and\n", + " return the fidelity and gradient of the fidelity to the global target_derivative\n", + " \"\"\"\n", + "\n", + " # Reshape flat parameter list to form accepted by state_gradient: [[hx0,hz0],[hx1,hz1,]...]\n", + " reshapedparas = [i for i in (paras.reshape((-1,num_params))).tolist() for j in range(2)]\n", + " reshapedparas = np.array(reshapedparas)\n", + "\n", + " gradient_dict = oqupy.state_gradient(\n", + " system=parametrized_system,\n", + " initial_state=initial_state,\n", + " target_derivative=target_derivative,\n", + " process_tensors=[process_tensor],\n", + " parameters=reshapedparas,\n", + " num_steps=num_steps,\n", + " progress_type='silent')\n", + " \n", + " fs=gradient_dict['final_state']\n", + " gps=gradient_dict['gradient']\n", + " fidelity=np.sum(fs*target_derivative.T)\n", + "\n", + " # Adding adjacent elements\n", + " for i in range(0,gps.shape[0],2): \n", + " gps[i,:]=gps[i,:]+gps[i+1,:]\n", + " \n", + " gps=gps[0::2]\n", + "\n", + " # Return the minus the gradient as infidelity is being minimized \n", + " return 1-fidelity.real,(-1.0*gps.reshape((-1)).real).tolist()\n", + "\n", + "# Set upper and lower bounds on control parameters\n", + "x_bound = [-5*np.pi,5*np.pi]\n", + "y_bound = [0,0] \n", + "z_bound = [-np.pi,np.pi]\n", + "\n", + "bounds = np.zeros((num_steps*num_params,2))\n", + "\n", + "for i in range(0, num_params*num_steps,num_params):\n", + " bounds[i] = x_bound\n", + " bounds[i+1] = y_bound\n", + " bounds[i+2] = z_bound\n", + "\n", + "y0 = np.zeros(num_steps)\n", + "z0 = np.ones(num_steps) * (np.pi) / (dt*num_steps)\n", + "x0 = np.zeros(num_steps)\n", + "\n", + "# Flatten list for input to optimizer\n", + "parameter_list=[item for pair in zip(x0, y0, z0) for item in pair]\n", + "\n", + "optimization_result = minimize(\n", + " fun=infidandgrad,\n", + " x0=parameter_list,\n", + " method='L-BFGS-B',\n", + " jac=True,\n", + " bounds=bounds,\n", + " options = {'disp':True, 'gtol': 7e-05}\n", + ")\n", + "\n", + "print(\"The maximal fidelity was found to be : \",1-optimization_result.fun)\n", + "\n", + "#print(\"The Jacobian was found to be : \",optimization_result.jac)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "\n", + "optimized_params = optimization_result.x\n", + "reshapedparas=np.array([i for i in (optimized_params.reshape((-1,num_params))).tolist() for j in range(2)])\n", + "\n", + "# Input optimized controls into state_gradient to show dynamics of system under optimized fields\n", + "optimized_dynamics = state_gradient(\n", + " system=parametrized_system,\n", + " initial_state=initial_state,\n", + " target_derivative=target_derivative,\n", + " process_tensors=[process_tensor],\n", + " parameters=reshapedparas,\n", + " num_steps=num_steps)\n", + "\n", + "fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1) \n", + "fig.suptitle(\"Optimisation results\")\n", + "\n", + "field_labels = [\"x\",\"y\",\"z\"]\n", + "for i in range(0,num_params):\n", + " ax1.plot(t[:-1],optimization_result['x'][i::num_params],label=field_labels[i])\n", + "ax1.set_ylabel(r\"$h_i$\",rotation=0,fontsize=16)\n", + "ax1.set_xlabel(\"t\")\n", + "ax1.legend()\n", + "\n", + "dynamics = optimized_dynamics['dynamics']\n", + "\n", + "t, bloch_x =optimized_dynamics['dynamics'].expectations(op.sigma(\"x\"))\n", + "t, bloch_y = optimized_dynamics['dynamics'].expectations(op.sigma(\"y\"))\n", + "t, bloch_z = optimized_dynamics['dynamics'].expectations(op.sigma(\"z\"))\n", + "\n", + "ax2.plot(t,bloch_x,label='x')\n", + "ax2.plot(t,bloch_y,label='y')\n", + "ax2.plot(t,bloch_z,label='z')\n", + "bloch_length = np.sqrt(bloch_x**2 +bloch_y**2 + bloch_z**2)\n", + "ax2.plot(t,bloch_length,label=r'$|\\mathbf{\\sigma}|$')\n", + "ax2.set_ylabel(r\"$\\langle \\sigma \\rangle$\",rotation=0,fontsize=16)\n", + "ax2.set_xlabel(\"t\")\n", + "\n", + "plt.legend()\n", + "plt.show()\n", + "\n", + "# -----------------------------------------------------------------------------\n", + "\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "oqupy", + "language": "python", + "name": "oqupy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.19" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/fidelity-gradient.py b/examples/fidelity-gradient.py index 32da49a0..32e7f33e 100644 --- a/examples/fidelity-gradient.py +++ b/examples/fidelity-gradient.py @@ -23,6 +23,10 @@ num_steps = 50 dt = end_time / num_steps +# this is the length of the PT to compute, which may be longer than the time above, to which the dynamics is considered + +pt_end_time = 5.0 + # -- bath -- alpha = 0.126 omega_cutoff = 3.04 @@ -59,7 +63,7 @@ process_tensor = oqupy.pt_tempo_compute( bath=bath, start_time=0.0, - end_time=end_time, + end_time=pt_end_time, parameters=pt_tempo_parameters) # --- Define parameterized system ---------------------------------------------- @@ -83,7 +87,8 @@ def hamiltonian(x, y, z): initial_state=initial_state, target_derivative=target_derivative, process_tensors=[process_tensor], - parameters=parameters) + parameters=parameters, + num_steps=num_steps) final_state = fidelity_dict['dynamics'].states[-1] v_final_state = np.reshape(final_state,hs_dim**2) @@ -126,6 +131,7 @@ def infidandgrad(paras): target_derivative=target_derivative, process_tensors=[process_tensor], parameters=reshapedparas, + num_steps=num_steps, progress_type='silent') fs=gradient_dict['final_state'] @@ -182,7 +188,8 @@ def infidandgrad(paras): initial_state=initial_state, target_derivative=target_derivative, process_tensors=[process_tensor], - parameters=reshapedparas) + parameters=reshapedparas, + num_steps=num_steps) fig, (ax1, ax2) = plt.subplots(nrows=2, ncols=1) fig.suptitle("Optimisation results") diff --git a/examples/tti-optimization-example.ipynb b/examples/tti-optimization-example.ipynb new file mode 100644 index 00000000..d312d87f --- /dev/null +++ b/examples/tti-optimization-example.ipynb @@ -0,0 +1,189 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.path.insert(0,'..')\n", + "\n", + "import oqupy\n", + "import oqupy.operators as op\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "import time\n", + "\n", + "import numpy as np\n", + "from scipy.integrate import solve_ivp\n", + "from scipy.interpolate import interp1d\n", + "from scipy.optimize import minimize,Bounds\n", + "\n", + "pt_parameters = {'tmax':0.5,\n", + " 'steps':50, \n", + " 'dkmax':50,\n", + " 'epsrel':10**(-5),\n", + " 'alpha':0.05,\n", + " 'omega_cutoff':2.0*np.pi*5.0, \n", + " 'temp':0.0}\n", + "\n", + "\n", + "opt_parameters = {'x0': pt_parameters['omega_cutoff'],\n", + " 'z0': 0.0,\n", + " 'lbx':2.0*np.pi*0.0,\n", + " 'ubx':2.0*np.pi*10.0,\n", + " 'lbz':0.0,\n", + " 'ubz':0.0,\n", + " 'initial_state':'mixed',\n", + " 'target_state':'x-'}\n", + "\n", + "\n", + "omega_cutoff = pt_parameters['omega_cutoff']\n", + "alpha = pt_parameters['alpha']\n", + "temperature = pt_parameters['temp']\n", + "epsrel = pt_parameters['epsrel']\n", + "t_max = pt_parameters['tmax']\n", + "# dt = 1./omega_cutoff/np.sqrt(3)\n", + "num_steps = pt_parameters['steps']\n", + "dt = t_max/num_steps\n", + "dkmax = pt_parameters['dkmax']\n", + "\n", + "initial_state = op.spin_dm('mixed')\n", + "target_state = op.spin_dm('x-')\n", + "target_derivative = target_state.T\n", + "\n", + "t_list = np.linspace(0,t_max,num_steps+1)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "import oqupy as oqupy\n", + "from oqupy.backends.itebd_tempo import iTEBD_TEMPO_oqupy\n", + "from oqupy.process_tensor import TTInvariantProcessTensor\n", + "from oqupy.tti_tempo import TTITempo\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt\n", + "\n", + "correlations = oqupy.PowerLawSD(alpha=alpha,\n", + " zeta=1,\n", + " cutoff=omega_cutoff,\n", + " cutoff_type='exponential',\n", + " temperature=temperature)\n", + "bath = oqupy.Bath(op.sigma(\"z\")/2.0, correlations)\n", + "parameters=oqupy.TempoParameters(dt=dt,epsrel=epsrel,dkmax=dkmax)\n", + "#pt=oqupy.pt_tempo_compute(\n", + "# bath=bath,\n", + "# start_time=0.0,\n", + "# end_time=t_max,\n", + "# parameters=parameters)\n", + "pt=TTITempo(bath,start_time=0.0,parameters=parameters)\n", + "process_tensor_tebd = pt.get_process_tensor()" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def hamiltonian_t(t):\n", + " return pt_parameters['omega_cutoff']*op.sigma('x')/2 + 0.0*op.sigma('z')/2\n", + "\n", + "system = oqupy.TimeDependentSystem(hamiltonian_t)\n", + "\n", + "Rho_0=oqupy.operators.spin_dm('mixed')\n", + "\n", + "dynamics = oqupy.compute_dynamics(\n", + " process_tensor=process_tensor_tebd, \n", + " system=system,\n", + " initial_state=Rho_0,\n", + " start_time=0,\n", + " num_steps=num_steps)\n", + "\n", + "t, s_x = dynamics.expectations(op.sigma('x'), real=True)\n", + "\n", + "plt.plot(t,np.log10(np.array(s_x)+1),label='sx')\n", + "plt.ylim((-5.0,1.0))\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "def discrete_hamiltonian(hx,hz):\n", + " return hx*op.sigma('x')/2 + hz*op.sigma('z')/2\n", + "system = oqupy.ParameterizedSystem(discrete_hamiltonian)\n", + "\n", + "h_z = np.ones(2*num_steps)*0 \n", + "h_x = np.ones(2*num_steps)*omega_cutoff\n", + "\n", + "parameters = np.vstack((h_x,h_z)).T\n", + "parameters.shape\n", + "\n", + "\n", + "grad_res = oqupy.state_gradient(\n", + " system=system,\n", + " initial_state=initial_state,\n", + " target_derivative=target_derivative,\n", + " process_tensors=[process_tensor_tebd],\n", + " num_steps=30,\n", + " parameters=parameters)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [ + "# some code for behind the scenes use to see what comes back from the gradient calculation\n", + "\n", + "grad_prop, dynamics = oqupy.compute_gradient_and_dynamics(\n", + " system=system,\n", + " initial_state=initial_state,\n", + " target_derivative=target_derivative,\n", + " process_tensors=[process_tensor_tebd],\n", + " parameters=parameters,\n", + " start_time=0.0,\n", + " dt=dt,\n", + " num_steps=27,\n", + " progress_type=None)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "oqupy", + "language": "python", + "name": "oqupy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.19" + } + }, + "nbformat": 4, + "nbformat_minor": 4 +} diff --git a/examples/tti-tempo.ipynb b/examples/tti-tempo.ipynb new file mode 100644 index 00000000..89e89801 --- /dev/null +++ b/examples/tti-tempo.ipynb @@ -0,0 +1,125 @@ +{ + "cells": [ + { + "cell_type": "code", + "execution_count": null, + "id": "8a85fd69", + "metadata": {}, + "outputs": [], + "source": [ + "# Demonstration of iTEBD within OQuPy\n", + "# Uses the original iTEBD code by Valentin Link from https://github.com/val-link/iTEBD-TEMPO.git.\n", + "# Minimally modified by Paul Eastham to take a bath correlation object from OQuPy.\n", + "# The resulting object can turned into a new proces tensor subclass, TTInvariantProcessTensor, which should function as all other PTs in OQuPy\n", + "# Todo: fully integrate the creation of iTEBD process tensors so we can skip calls to iTEBD_TEMPO_oqupy class." + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f906bed1", + "metadata": {}, + "outputs": [], + "source": [ + "import sys\n", + "sys.path.insert(0,'..')\n", + "\n", + "import oqupy as oqupy\n", + "from oqupy.backends.itebd_tempo import iTEBD_TEMPO_oqupy\n", + "from oqupy.process_tensor import TTInvariantProcessTensor\n", + "from oqupy.tti_tempo import TTITempo\n", + "import numpy as np\n", + "import matplotlib.pyplot as plt" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "f314d374", + "metadata": {}, + "outputs": [], + "source": [ + "sigma_z = oqupy.operators.sigma(\"z\")\n", + "sigma_x = oqupy.operators.sigma(\"x\")\n", + "sigma_y = oqupy.operators.sigma(\"y\")\n", + "omega_cutoff = 3.04 \n", + "alpha = 0.126\n", + "temperature = 0.1309\n", + "correlations = oqupy.PowerLawSD(alpha=alpha,\n", + " zeta=3,\n", + " cutoff=omega_cutoff,\n", + " cutoff_type='gaussian',\n", + " temperature=temperature)\n", + "bath = oqupy.Bath(sigma_x/2.0, correlations)\n", + "parameters=oqupy.TempoParameters(dt=0.1,epsrel=1e-7,dkmax=100)\n", + "pt=TTITempo(bath,start_time=0.0,parameters=parameters)\n" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "6dc20eab", + "metadata": {}, + "outputs": [], + "source": [ + "detuning = lambda t: 0.0 * t\n", + "def gaussian_shape(t, area = 1.0, tau = 1.0, t_0 = 0.0):\n", + " return area/(tau*np.sqrt(np.pi)) * np.exp(-(t-t_0)**2/(tau**2))\n", + "\n", + "def hamiltonian_t(t):\n", + " return detuning(t)/2.0 * sigma_z \\\n", + " + gaussian_shape(t, area = np.pi/2.0, tau = 0.5)/2.0 * sigma_y\n", + "\n", + "system = oqupy.TimeDependentSystem(hamiltonian_t)\n", + "\n", + "Rho_0=oqupy.operators.spin_dm('x-')\n", + "\n", + "dynamics = oqupy.compute_dynamics(\n", + " process_tensor=pt.get_process_tensor(), \n", + " system=system,\n", + " initial_state=Rho_0,\n", + " start_time=-5.0,\n", + " num_steps=200)\n", + "\n", + "t, s_x = dynamics.expectations(sigma_x, real=True)\n", + "_, s_y = dynamics.expectations(sigma_y, real=True)\n", + "_, s_z = dynamics.expectations(sigma_z, real=True)\n", + "\n", + "plt.plot(t,s_x,label='sx')\n", + "plt.plot(t,s_y,label='sy')\n", + "plt.plot(t,s_z,label='sz')\n", + "plt.ylim((-1.0,1.0))\n", + "plt.legend(loc=4)" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "be7831c4-0a8a-4625-b9d9-da735ed1210d", + "metadata": {}, + "outputs": [], + "source": [] + } + ], + "metadata": { + "kernelspec": { + "display_name": "oqupy", + "language": "python", + "name": "oqupy" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.10.19" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/oqupy/__init__.py b/oqupy/__init__.py index 1a72fe52..2c4e99bf 100644 --- a/oqupy/__init__.py +++ b/oqupy/__init__.py @@ -39,13 +39,14 @@ (https://doi.org/10.1103/PhysRevLett.129.173001) (2022). [8] Fux et al., [Phys. Rev. Lett. 126, 200401] (https://doi.org/10.1103/PhysRevLett.126.200401) (2021). -[9] Butler et al., [Phys. Rev. Lett. 132, 060401 ] +[9] Butler et al., [Phys. Rev. Lett. 132, 060401] (https://doi.org/10.1103/PhysRevLett.132.060401}) (2024). [10] Gribben et al., [Quantum, 6, 847] (https://doi.org/10.22331/q-2022-10-25-847) (2022). [11] Chiu et al., [Phys. Rev. A 106, 012204] - (https://doi.org/10.1103/PhysRevA.106.012204}) (2022). - + (https://doi.org/10.1103/PhysRevA.106.012204) (2022). +[12] Link et al., [Phys. Rev. Lett. 132, 200403] + (https://doi.org/10.1103/PhysRevLett.132.200403) (2024). """ from oqupy.version import __version__ @@ -80,7 +81,11 @@ 'PtTebdParameters', 'PtTempo', 'pt_tempo_compute', + 'TTITempo', + 'tti_tempo_compute', 'SimpleProcessTensor', + 'SimpleProcessTensorFinite', + 'SimpleProcessTensorInfinite', 'state_gradient', 'System', 'SystemChain', @@ -126,6 +131,8 @@ from oqupy.process_tensor import import_process_tensor from oqupy.process_tensor import TrivialProcessTensor from oqupy.process_tensor import SimpleProcessTensor +from oqupy.process_tensor import SimpleProcessTensorFinite +from oqupy.process_tensor import SimpleProcessTensorInfinite from oqupy.process_tensor import FileProcessTensor from oqupy.pt_tebd import PtTebd @@ -141,6 +148,9 @@ from oqupy.pt_tempo import PtTempo from oqupy.pt_tempo import pt_tempo_compute +from oqupy.tti_tempo import TTITempo +from oqupy.tti_tempo import tti_tempo_compute + from oqupy.tempo import Tempo from oqupy.tempo import TempoParameters from oqupy.tempo import GibbsTempo diff --git a/oqupy/backends/tti_tempo_backend.py b/oqupy/backends/tti_tempo_backend.py new file mode 100644 index 00000000..9e778273 --- /dev/null +++ b/oqupy/backends/tti_tempo_backend.py @@ -0,0 +1,208 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Module for the time translational invariant time evolving matrix product +operator algorithm (TTI-TEMPO) backend. This module is based on [Link2024]. + +**[Link2024]** +V. Link, H. Tu, and W. T. Strunz, *Open Quantum System Dynamics from Infinite +Tensor Network Contraction*, `Phys. Rev. Lett. 132, 200403 +`__ (2024). + +Original code taken from https://github.com/val-link/iTEBD-TEMPO.git +# Implementation of iTEBD-TEMPO +# Author: Valentin Link (valentin.link@tu-dresden.de) +# Original code from https://github.com/val-link/iTEBD-TEMPO.git +# Modified by Paul Eastham (easthamp@tcd.ie) to use the OQuPy BathCorrelations +# class to define the bath correlations rather than the bath correlation +# function itself. +# Please cite the corresponding publication: +# https://doi.org/10.1103/PhysRevLett.132.200403. +""" + +from typing import Callable, Dict + +import numpy as np +from scipy.linalg import norm, svd +from scipy.sparse.linalg import eigs + +from oqupy.process_tensor import SimpleProcessTensorInfinite + +class TTITempoBackend(): + """ + Base class for TTI process tensor tempo backends. + + Parameters + ---------- + influence: callable(int) -> ndarray + Callable that takes an integer `step` and returns the influence super + operator of that `step`. + process_tensor: BaseProcessTensor + Todo + dkmax: int + Number of influences to include. If ``dkmax == None`` then all + influences are included. + epsrel: float + Maximal relative SVD truncation error. + rank: int + Maximal SVD truncation rank + config: Dict + Todo + """ + + def __init__( + self, + dimension: int, + influence: Callable[[int], np.ndarray], + process_tensor: SimpleProcessTensorInfinite, + dkmax: int, + epsrel: float, + rank: int, + config: Dict): + + self.repeating_cell = (1, 1) + self._dkmax = dkmax + self._dimension = dimension + self._influence = influence + self._process_tensor = process_tensor + self._nu_dim = self._dimension**2 + 1 + self._epsrel = epsrel + self._rank = rank + self._itol = config["itol"] + self._A = None + self._B = None + self._sAB = None + self._sBA = None + self._rank_is_one = None + self._step = None + + @property + def step(self) -> int: + """The current step in the TTI-TEMPO computation. """ + return self._step + + @property + def num_steps(self) -> int: + """The current step in the TTI-TEMPO computation. """ + return self._dkmax + + def initialize(self) -> None: + """Initializes the iTEBD algorithm. """ + self._A = np.ones((1, self._nu_dim, 1)) + self._B = np.ones((1, self._nu_dim, 1)) + self._sAB = np.ones((1)) + self._sBA = np.ones((1)) + self._rank_is_one = True + self._step = 1 + + def compute_step(self) -> None: + """ + Compute a step of the iTEBD algorithm. + """ + i_tens = self._influence(self._dkmax - self._step) + + if self._step % 2 == 0: + self._sBA, self._sAB = self._sAB, self._sBA + self._B, self._A, = self._A, self._B + + # renormalize weights + self._sAB = self._sAB * norm(self._sBA) + self._sBA = self._sBA / norm(self._sBA) + + # ensure weights are above tolerance (needed for inversion) + self._sBA[np.abs(self._sBA) < self._itol] = self._itol + + # MPS - gate contraction + d1 = i_tens.shape[1] + d2 = i_tens.shape[-1] + rank_BA = self._sBA.shape[0] + + if self._step == self._dkmax: + d1 = 1 + u, s_vals, v = svd(np.einsum('a,acd,d,dcg,g,i,c->aicg', self._sBA, + self._A, self._sAB, self._B, self._sBA, + np.ones((1)), np.diagonal(i_tens) + ).reshape([d1*rank_BA, d2*rank_BA]), + full_matrices=False) + else: + u, s_vals, v = svd(np.einsum('a,acd,d,dfg,g,fc->afcg', self._sBA, + self._A, self._sAB, self._B, self._sBA, + i_tens + ).reshape([d1*rank_BA, d2*rank_BA]), + full_matrices=False) + + # truncate singular values + if self._epsrel is None: + rank_new = min(self._rank, len(s_vals)) + else: + s_vals_sum = np.cumsum(s_vals) / np.sum(s_vals) + rank_rtol = np.searchsorted(s_vals_sum, 1 - self._epsrel) + 1 + rank_new = min(self._rank, len(s_vals), rank_rtol) + u = u[:, :rank_new].reshape(self._sBA.shape[0], d1 * rank_new) + v = v[:rank_new, :].reshape(rank_new * d2, rank_BA) + + # factor out sAB weights from A and B + self._A = (np.diag(1 / self._sBA) @ u).reshape( + self._sBA.shape[0], d1, rank_new) + self._B = (v @ np.diag(1 / self._sBA)).reshape(rank_new, d2, rank_BA) + + # new weights + self._sAB = s_vals[:rank_new] + + if self._step % 2 == 0: + self._sBA, self._sAB = self._sAB, self._sBA + self._B, self._A, = self._A, self._B + + if self._rank_is_one: + if all([self._sAB.shape[0] == 1, self._sAB.shape[-1] == 1, + self._sBA.shape[0] == 1, self._sBA.shape[-1] == 1]): + # reset to initial mps if rank is still one + self._sAB = np.ones((1)) + self._sBA = np.ones((1)) + self._A = np.ones((1, self._nu_dim, 1)) + self._B = np.ones((1, self._nu_dim, 1)) + else: + self._rank_is_one = False +# print(f"Effective memory depth dkmax={self._dkmax - self._step + 1}") + if self._step == 1: + print("Warning: the memory cutoff dkmax may be too small"\ + "for the given epsrel value. The algorithm may"\ + "become unstable and inaccurate. It is recommended"\ + "to increase dkmax until this message does no"\ + "longer appear.") + self._step += 1 + return self._step-1 < self.num_steps + + def update_process_tensor(self) -> None: + """Update the process tensor. """ + assert self._step >= self._dkmax + f = np.squeeze(np.einsum('i,ikl,l,lop->ikop', self._sAB, self._B, + self._sBA, self._A)) + + if self._sAB.shape[0] == 1: + # handle trivial f + f = np.ones((1, self._nu_dim, 1)) + + # compute f[:,-1,:]^\inf = v_r * v_l^T using Lanczos + _, v_r = eigs(f[:, -1, :], 1, which='LR') + _, v_l = eigs(f[:, -1, :].T, 1, which='LR') + v_l = v_l / (v_l[:, 0] @ v_r[:, 0]) + + self._process_tensor.set_mpo_tensor(0, np.einsum('ia,ikj->ajk', + v_l, f[:,:-1,:])) + self._process_tensor.set_mpo_tensor(1, np.swapaxes(f[:,:-1,:], 1, 2)) + self._process_tensor.set_cap_tensor(0, np.array([1.0])) + self._process_tensor.set_cap_tensor(1, v_r[:, 0]) + try: + self._process_tensor.repeating_cell = self.repeating_cell + except Exception: + pass diff --git a/oqupy/config.py b/oqupy/config.py index 592dfa7b..398528f9 100644 --- a/oqupy/config.py +++ b/oqupy/config.py @@ -68,6 +68,12 @@ PT_TEBD_DEFAULT_EPSREL = 1.0e-5 +# -- TTI_TEMPO --------------------------------------------------------------- + +# Default TTI-TEMPO backend configuration +TTI_TEMPO_BACKEND_CONFIG = {"itol": 1e-13} + + # -- BATH -------------------------------------------------------------------- diff --git a/oqupy/gradient.py b/oqupy/gradient.py index 9c177436..a4aab3af 100644 --- a/oqupy/gradient.py +++ b/oqupy/gradient.py @@ -35,7 +35,9 @@ def state_gradient( process_tensors: List[BaseProcessTensor], parameters: ndarray, start_time: Optional[float] = 0.0, - progress_type: Optional[Text] = None) -> Dict: + num_steps: Optional[int]=None, + progress_type: Optional[Text] = None, + only_dynamics: Optional[bool] = False) -> Dict: """ Compute system dynamics and gradient of an objective function Z with respect to a parameterized System for a given set of control @@ -75,7 +77,6 @@ def state_gradient( """ check_isinstance(parameters, ndarray, 'parameters') - num_steps = len(process_tensors[0]) dt = process_tensors[0].dt num_parameters = parameters.shape[1] @@ -88,26 +89,34 @@ def state_gradient( start_time=start_time, dt=dt, num_steps=num_steps, - progress_type=progress_type) - - - get_half_props= system.get_propagators(dt,parameters) - get_prop_derivatives = system.get_propagator_derivatives(dt,parameters) - - final_derivs = _chain_rule( - adjoint_tensor=grad_prop, - dprop_dparam=get_prop_derivatives, - propagators=get_half_props, - num_steps=num_steps, - num_parameters=num_parameters, - progress_type=progress_type) - - return_dict = { - 'final_state':dynamics.states[-1], - 'gradprop':grad_prop, - 'gradient':final_derivs, - 'dynamics':dynamics - } + progress_type=progress_type, + only_dynamics=only_dynamics) + + if only_dynamics: + return_dict = { + 'final_state':dynamics.states[-1], + 'gradprop':None, + 'gradient':None, + 'dynamics':dynamics + } + else: + get_half_props= system.get_propagators(dt,parameters) + get_prop_derivatives = system.get_propagator_derivatives(dt,parameters) + + final_derivs = _chain_rule( + adjoint_tensor=grad_prop, + dprop_dparam=get_prop_derivatives, + propagators=get_half_props, + num_steps=len(grad_prop), + num_parameters=num_parameters, + progress_type=progress_type) + + return_dict = { + 'final_state':dynamics.states[-1], + 'gradprop':grad_prop, + 'gradient':final_derivs, + 'dynamics':dynamics + } return return_dict @@ -177,7 +186,8 @@ def compute_gradient_and_dynamics( num_steps: Optional[int] = None, control: Optional[Control] = None, record_all: Optional[bool] = True, - progress_type: Optional[Text] = None) -> Tuple[List, Dynamics]: + progress_type: Optional[Text] = None, + only_dynamics: Optional[bool] = False) -> Tuple[List, Dynamics]: """ Compute some objective function and calculate its gradient w.r.t. some control parameters, accounting for interaction with an environment @@ -211,7 +221,10 @@ def compute_gradient_and_dynamics( progress_type: str (default = None) The progress report type during the computation. Types are: {``silent``, ``simple``, ``bar``}. If `None` then - the default progress type is used. + the default progress type is used.# + only_dynamics: bool (default = False) + Set to true to compute only the dynamics, target_derivative is not + used in this case. Returns: -------- @@ -242,6 +255,9 @@ def compute_gradient_and_dynamics( num_envs = len(process_tensors) + if num_steps is None: + num_steps=len(process_tensors[0]) + # -- prepare propagators -- propagators = system.get_propagators(dt, parameters) @@ -261,6 +277,7 @@ def controls(step: int): # Initial state including the bond legs to the environments with: # edges 0, 1, .., num_envs-1 are the bond legs of the environments # edge -1 is the state leg + initial_ndarray = initial_state.reshape(hs_dim**2) initial_ndarray.shape = tuple([1]*num_envs+[hs_dim**2]) current_node = tn.Node(initial_ndarray) @@ -333,6 +350,9 @@ def controls(step: int): dynamics = Dynamics(times=list(times),states=states) + if only_dynamics: + return None,dynamics + # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ # ~~~~~ Backpropagation ~~~~~~ # ~~~~~~~~~~~~~~~~~~~~~~~~~~~~ @@ -353,7 +373,8 @@ def controls(step: int): target_ndarray = target_derivative target_ndarray = target_ndarray.reshape(hs_dim**2) - target_ndarray.shape = tuple([1]*num_envs+[hs_dim**2]) + # target_ndarray.shape = tuple([1]*num_envs+[hs_dim**2]) + target_ndarray = np.outer(caps,target_ndarray) current_node = tn.Node(target_ndarray) current_edges = current_node[:] diff --git a/oqupy/process_tensor.py b/oqupy/process_tensor.py index d5df0623..14d8af3a 100644 --- a/oqupy/process_tensor.py +++ b/oqupy/process_tensor.py @@ -35,6 +35,9 @@ from oqupy import util from oqupy.version import __version__ +NoneType = type(None) + + class BaseProcessTensor(BaseAPIClass, ABC): """ Abstract base class for process tensors in matrix product operator form @@ -275,6 +278,26 @@ def __len__(self) -> int: """Length of process tensor. """ return len(self._mpo_tensors) + @property + def repeating_cell(self) -> Union[tuple[int, int], NoneType]: + """ + Repeating unit cell of a time translational invariant process tensor. + """ + return None + + def effective_step(self, + step: int) -> None: + """ + Get the effective step for the _mpo_tensors and _cap_tensors lists + when there's a repeating cell. + """ + if self.repeating_cell is None: + k = step + else: + n, l = self.repeating_cell[0], self.repeating_cell[1] + k = step if step < n else n + ((step-n) % l) + return k + def set_initial_tensor( self, initial_tensor: Optional[ndarray] = None) -> None: @@ -339,12 +362,13 @@ def get_mpo_tensor( Applies the transformation (stored in `.transform_in` and `.transform_out`) when `transformed` is true. """ + step = self.effective_step(step) length = len(self._mpo_tensors) if step >= length or step < 0: raise IndexError("Process tensor index out of bound. ") tensor = self._mpo_tensors[step] if len(tensor.shape) == 3: - tensor = util.create_delta(tensor, [0, 1, 2, 2]) + tensor = util.create_delta_lastindex(tensor) if transformed is False: return tensor if self._transform_in is not None: @@ -359,6 +383,7 @@ def get_cap_tensor(self, step: int) -> ndarray: """ Get the cap tensor (vector) to terminate the PT-MPO at time step `step`. """ + step = self.effective_step(step) length = len(self._cap_tensors) if step >= length or step < 0: return None @@ -377,6 +402,37 @@ def get_bond_dimensions(self) -> ndarray: bond_dims.append(self._mpo_tensors[-1].shape[1]) return np.array(bond_dims) + def export(self, filename: Text, overwrite: bool = False): + """Export the process tensor as a FileProcessTensor.""" + if overwrite: + mode = "overwrite" + else: + mode = "write" + + pt_file = FileProcessTensor( + mode=mode, + filename=filename, + hilbert_space_dimension=self._hs_dim, + dt=self._dt, + transform_in=self._transform_in, + transform_out=self._transform_out, + name=self.name, + description=self.description) + + pt_file.set_initial_tensor(self._initial_tensor) + for step, mpo in enumerate(self._mpo_tensors): + pt_file.set_mpo_tensor(step, mpo) + for step, cap in enumerate(self._cap_tensors): + pt_file.set_cap_tensor(step, cap) + pt_file.close() + +class SimpleProcessTensorFinite(SimpleProcessTensor): + """Finite Process Tensor""" + + @property + def repeating_cell(self) -> NoneType: + return None + def compute_caps(self) -> None: """ Compute and store all caps from the PT-MPO. @@ -405,29 +461,22 @@ def compute_caps(self) -> None: last_cap = new_cap self._cap_tensors = caps - def export(self, filename: Text, overwrite: bool = False): - """Export the process tensor as a FileProcessTensor.""" - if overwrite: - mode = "overwrite" - else: - mode = "write" - pt_file = FileProcessTensor( - mode=mode, - filename=filename, - hilbert_space_dimension=self._hs_dim, - dt=self._dt, - transform_in=self._transform_in, - transform_out=self._transform_out, - name=self.name, - description=self.description) +class SimpleProcessTensorInfinite(SimpleProcessTensor): + """ + Simple infinite process tensor in matrix product operator form with + a repeating unit cell. + """ - pt_file.set_initial_tensor(self._initial_tensor) - for step, mpo in enumerate(self._mpo_tensors): - pt_file.set_mpo_tensor(step, mpo) - for step, cap in enumerate(self._cap_tensors): - pt_file.set_cap_tensor(step, cap) - pt_file.close() + @property + def max_step(self) -> Union[int, float]: + """Maximal number of time steps.""" + return float('inf') + + @property + def repeating_cell(self) -> tuple[int, int]: + """Repeating unit cell of an infinite process tensor.""" + return (1, 1) HDF5None = [np.nan] @@ -836,19 +885,23 @@ def import_process_tensor( process_tensor_type: Text Type of process tensor object to create. May be 'file' to create `FileProcessTensor`, - or 'simple', to create a `SimpleProcessTensor`. + or 'simple', to create a generic `SimpleProcessTensor`, + or 'simple-finite', to create a `SimpleProcessTensorFinite`, + or 'simple-infinite' to create a `SimpleProcessTensorInfinite`. Returns ------- process_tensor: BaseProcessTensor The process tensor object with the data from the .hdf5 file. """ + map_dict = {"simple": SimpleProcessTensor, + "simple-finite": SimpleProcessTensorFinite, + "simple-infinite": SimpleProcessTensorInfinite} pt_file = FileProcessTensor(mode="read", filename=filename) - if process_tensor_type is None or process_tensor_type == "file": pt = pt_file - elif process_tensor_type == "simple": - pt = SimpleProcessTensor( + elif process_tensor_type in map_dict: + pt = map_dict[process_tensor_type]( hilbert_space_dimension=pt_file.hilbert_space_dimension, dt=pt_file.dt, transform_in=pt_file.transform_in, diff --git a/oqupy/pt_tempo.py b/oqupy/pt_tempo.py index 9cc6f822..6e3c6d06 100644 --- a/oqupy/pt_tempo.py +++ b/oqupy/pt_tempo.py @@ -47,7 +47,7 @@ from oqupy.config import PT_DEFAULT_TOLERANCE from oqupy.config import PT_TEMPO_BACKEND_CONFIG from oqupy.process_tensor import BaseProcessTensor -from oqupy.process_tensor import SimpleProcessTensor +from oqupy.process_tensor import SimpleProcessTensorFinite from oqupy.process_tensor import FileProcessTensor from oqupy.backends.pt_tempo_backend import PtTempoBackend from oqupy.tempo import TempoParameters @@ -57,7 +57,7 @@ from oqupy.util import get_progress -PT_CLASS = {"simple": SimpleProcessTensor} +PT_CLASS = {"simple": SimpleProcessTensorFinite} @@ -165,7 +165,7 @@ def _init_simple_process_tensor(self): else: transform_in = None transform_out = None - self._process_tensor = SimpleProcessTensor( + self._process_tensor = SimpleProcessTensorFinite( hilbert_space_dimension=self._dimension, dt=self._parameters.dt, transform_in=transform_in, diff --git a/oqupy/system_dynamics.py b/oqupy/system_dynamics.py index 1ddc9659..eb314e75 100644 --- a/oqupy/system_dynamics.py +++ b/oqupy/system_dynamics.py @@ -521,7 +521,8 @@ def _compute_dynamics_input_parse( "One of the elements in `process_tensor` is not of type " \ + "`{BaseProcessTensor.__name__}`.") if pt.get_initial_tensor() is not None: - raise NotImplementedError() + warnings.warn("Initial correlated states are not yet "\ + "implemented, ignoring the initial correlations.") check_true( hs_dim == pt.hilbert_space_dimension, "All process tensor must have the same Hilbert "\ @@ -795,6 +796,7 @@ def compute_correlations_nt( ops_times: List[Union[Indices, float, Tuple[float, float]]], ops_order: List[Text], initial_state: Optional[ndarray] = None, + num_steps: Optional[int] = None, start_time: Optional[float] = 0.0, dt: Optional[float] = None, progress_type: Text = None, @@ -828,6 +830,8 @@ def compute_correlations_nt( Initial system state. start_time: float (default = 0.0) Initial time. + num_steps: int + Optional number of time steps to be computed. dt: float (default = None) Time step size. progress_type: str (default = None) @@ -875,8 +879,16 @@ def compute_correlations_nt( #Input parsing; ensures that specified times that are not an integer multiple #of dt are assigned the closest integer multiple.------------------------------ + max_step = process_tensor.max_step - max_step = len(process_tensor) + if num_steps is None: + max_step = process_tensor.max_step + check_true( + max_step < np.inf, + "Variable `num_steps` must be specified because all "\ + "process tensors involved are infinite.") + else: + max_step = num_steps ops_times_=[] ret_times=[] #These are the times returned by the function @@ -1102,7 +1114,16 @@ def _parse_times(times, max_step, dt, start_time): ret_times = np.array([times]) elif isinstance(times, (slice, list)): try: - ret_times = np.arange(max_step + 1)[times] + if isinstance(times, slice): + max_time = max(i for i in [times.stop, times.start] + if i is not None) + ret_times = range(max_time+1)[times] + else: + ret_times = times + ret_times = np.fromiter(ret_times, dtype=np.int64, + count=len(ret_times)) + if ((ret_times < 0) | (ret_times > max_step)).any(): + raise IndexError except Exception as e: raise IndexError("Specified times are invalid or out of bound.") \ from e @@ -1122,8 +1143,11 @@ def _parse_times(times, max_step, dt, start_time): if index_end < 0 or index_end > max_step: raise IndexError("Specified end time is out of bound.") direction = 1 if index_start <= index_end else -1 - ret_times = np.arange( - max_step + 1)[index_start:index_end+direction:direction] + trange = range(max(index_start, index_end)+1) + ret_times = trange[index_start:index_end+direction:direction] + ret_times = np.fromiter(ret_times, dtype=np.int64, count=len(ret_times)) + if ((ret_times < 0) | (ret_times > max_step)).any(): + raise IndexError else: raise TypeError("Parameters `times_a` and `times_b` must be either " \ + "int, slice, list, or tuple.") diff --git a/oqupy/tempo.py b/oqupy/tempo.py index 896679bc..d833f28d 100644 --- a/oqupy/tempo.py +++ b/oqupy/tempo.py @@ -47,6 +47,8 @@ from oqupy.util import check_convert, check_isinstance, check_true,\ get_progress +NoneType = type(None) + class TempoParameters(BaseAPIClass): r""" Parameters for the TEMPO computation. @@ -63,6 +65,10 @@ class TempoParameters(BaseAPIClass): in the underlying tensor network algorithm). - It must be small enough such that the numerical compression (using tensor network algorithms) does not truncate relevant correlations. + rank: int (default = np.inf) + The maximal rank in the singular value truncation (done in the + underlying tensor network algorithm). The relative error in the + singular value truncation might end up higher than `epsrel`. tcut: float (default = None) Length of time :math:`t_\mathrm{cut}` included in the non-Markovian memory. - This should be large enough to capture all non-Markovian @@ -94,6 +100,7 @@ def __init__( self, dt: float, epsrel: float, + rank: Optional[Union[int, float, NoneType]] = np.inf, tcut: Optional[float] = None, dkmax: Optional[int] = None, add_correlation_time: Optional[float] = None, @@ -119,6 +126,17 @@ def __init__( raise ValueError("Argument 'epsrel' must be positive.") self._epsrel = tmp_epsrel + try: + if rank == np.inf or rank is None: + tmp_rank = np.inf + else: + tmp_rank = int(rank) + except Exception as e: + raise TypeError("Argument 'rank' must be an integer.") from e + if tmp_rank <= 0: + raise ValueError("Argument 'rank' must be strictly positive.") + self._rank = tmp_rank + self._tcut, self._dkmax = _parameter_memory_input_parse( tcut, dkmax, dt) @@ -167,6 +185,7 @@ def __str__(self) -> Text: ret.append(" tcut [dkmax] = {} [{}] \n".format( self.tcut, self.dkmax)) ret.append(" epsrel = {} \n".format(self.epsrel)) + ret.append(" rank = {} \n".format(self.rank)) ret.append(" add_correlation_time = {} \n".format( self.add_correlation_time)) return "".join(ret) @@ -183,6 +202,11 @@ def epsrel(self) -> float: # epsrel is error tolerance for both correlation omega integration and svds + @property + def rank(self) -> float: + """The maximal rank in the singular value truncation.""" + return self._rank + @property def tcut(self) -> float: """Length of non-Markovian memory""" @@ -999,8 +1023,10 @@ def influence_matrix( delta=dt, time_1=time_1, time_2=time_2, - shape=shape, - epsrel=parameters.epsrel) + shape=shape) + # epsrel=parameters.epsrel) +#/!\ SVD truncation epsrel doesn't necessarily match the integration epsrel + op_p = coupling_acomm op_m = coupling_comm diff --git a/oqupy/tti_tempo.py b/oqupy/tti_tempo.py new file mode 100644 index 00000000..0bcd40fa --- /dev/null +++ b/oqupy/tti_tempo.py @@ -0,0 +1,264 @@ +# Licensed under the Apache License, Version 2.0 (the "License"); +# you may not use this file except in compliance with the License. +# You may obtain a copy of the License at +# +# http://www.apache.org/licenses/LICENSE-2.0 +# +# Unless required by applicable law or agreed to in writing, software +# distributed under the License is distributed on an "AS IS" BASIS, +# WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied. +# See the License for the specific language governing permissions and +# limitations under the License. +""" +Module for the time translational invariant time evolving matrix product +operator algorithm (TTI-TEMPO). This module is based on [Link2024]. + +**[Link2024]** +V. Link, H. Tu, and W. T. Strunz, *Open Quantum System Dynamics from Infinite +Tensor Network Contraction*, `Phys. Rev. Lett. 132, 200403 +`__ (2024). +""" + +from typing import Dict, Optional, Text, Union + +import numpy as np + +from oqupy.base_api import BaseAPIClass +from oqupy.bath import Bath +from oqupy.config import TTI_TEMPO_BACKEND_CONFIG +from oqupy.process_tensor import BaseProcessTensor, SimpleProcessTensorInfinite +from oqupy.process_tensor import FileProcessTensor +from oqupy.backends.tti_tempo_backend import TTITempoBackend +from oqupy.tempo import TempoParameters +from oqupy.tempo import influence_matrix +from oqupy.operators import left_right_super +from oqupy.util import get_progress + + +class TTITempo(BaseAPIClass): + """ + Class to facilitate a TTI-TEMPO computations. + + Parameters + ---------- + bath: Bath + The Bath (includes the coupling operator to the system). + parameters: TempoParameters + The parameters for the PT-TEMPO computation. + start_time: float + The start time. + backend_config: dict (default = None) + The configuration of the backend. If `backend_config` is + ``None`` then the default backend configuration is used. + name: str (default = None) + An optional name for the tempo object. + description: str (default = None) + An optional description of the tempo object. + """ + def __init__( + self, + bath: Bath, + start_time: float, + parameters: TempoParameters, + process_tensor_file: Optional[Union[Text, bool]] = None, + overwrite: Optional[bool] = False, + backend_config: Optional[Dict] = None, + name: Optional[Text] = None, + description: Optional[Text] = None) -> None: + """Create a TTITempo object. """ + assert isinstance(bath, Bath), \ + "Argument 'bath' must be an instance of Bath." + self._bath = bath + self._dimension = self._bath.dimension + self._correlations = self._bath.correlations + + super().__init__(name, description) + + try: + tmp_start_time = float(start_time) + except Exception as e: + raise AssertionError("Start time must be a float.") from e + self._start_time = tmp_start_time + + assert isinstance(parameters, TempoParameters), \ + "Argument 'parameters' must be an instance of TempoParameters." + self._parameters = parameters + + self._process_tensor = None + if process_tensor_file or isinstance(process_tensor_file, Text): + if isinstance(process_tensor_file, Text): + filename = process_tensor_file + else: + filename = None + self._init_file_process_tensor(filename, overwrite) + else: + self._init_infinite_process_tensor() + + if backend_config is None: + self._backend_config = TTI_TEMPO_BACKEND_CONFIG + else: + self._backend_config = TTI_TEMPO_BACKEND_CONFIG | backend_config + + self._coupling_comm = np.pad(self._bath._coupling_comm, [(0, 1)]) + self._coupling_acomm = np.pad(self._bath._coupling_acomm, [(0, 1)]) + + self._backend_instance = None + self._init_tti_tempo_backend() + + def _init_infinite_process_tensor(self): + """ToDo. """ + unitary = self._bath.unitary_transform + if not np.allclose(unitary, np.identity(self._dimension)): + transform_in = left_right_super(unitary.conjugate().T, + unitary).T + transform_out = left_right_super(unitary, + unitary.conjugate().T).T + else: + transform_in = None + transform_out = None + + self._process_tensor = SimpleProcessTensorInfinite( + hilbert_space_dimension=self._dimension, + dt=self._parameters.dt, + transform_in=transform_in, + transform_out=transform_out, + name=self.name, + description=self.description) + + def _init_file_process_tensor(self, filename, overwrite): + """ToDo. """ + unitary = self._bath.unitary_transform + if not np.allclose(unitary, np.identity(self._dimension)): + transform_in = left_right_super(unitary.conjugate().T, + unitary).T + transform_out = left_right_super(unitary, + unitary.conjugate().T).T + else: + transform_in = None + transform_out = None + + if overwrite: + mode = "overwrite" + else: + mode = "write" + self._process_tensor = FileProcessTensor( + mode=mode, + filename=filename, + hilbert_space_dimension=self._dimension, + dt=self._parameters.dt, + transform_in=transform_in, + transform_out=transform_out, + name=self.name, + description=self.description) + + def _init_tti_tempo_backend(self): + """Create and initialize the tti-tempo backend.""" + self._backend_instance = TTITempoBackend( + dimension=self._dimension, + influence=self._influence, + process_tensor=self._process_tensor, + dkmax=self._parameters.dkmax, + epsrel=self._parameters.epsrel, + rank=self._parameters.rank, + config=self._backend_config) + + def _influence(self, dk): + return influence_matrix( + dk, + parameters=self._parameters, + correlations=self._correlations, + coupling_acomm=self._coupling_acomm, + coupling_comm=self._coupling_comm) + + def compute(self, progress_type: Optional[Text] = None) -> None: + """ + Propagate (or continue to propagate) the TEMPO tensor network to + time `end_time`. + + Parameters + ---------- + progress_type: str (default = None) + The progress report type during the computation. Types are: + {``silent``, ``simple``, ``bar``}. If `None` then + the default progress type is used. + """ + if self._backend_instance.step is None: + self._backend_instance.initialize() + + progress = get_progress(progress_type) + title = "--> TTI-TEMPO computation:" + with progress(self._backend_instance.num_steps, title) as prog_bar: + while self._backend_instance.compute_step(): + prog_bar.update(self._backend_instance.step) + + def get_process_tensor( + self, + progress_type: Optional[Text] = None) -> BaseProcessTensor: + """ + Returns a the computed process tensor. It performs the computation if + it hasn't been already done. + + Parameters + ---------- + progress_type: str (default = None) + The progress report type during the computation. Types are: + {``silent``, ``simple``, ``bar``}. If `None` then + the default progress type is used. + + Returns + ------- + process_tensor: SimpleProcessTensorInfinite + The computed process tensor. + """ + if self._backend_instance.step is None or \ + self._backend_instance.step < self._backend_instance.num_steps: + self.compute(progress_type=progress_type) + + if len(self._process_tensor) < sum( + self._backend_instance.repeating_cell): + self._backend_instance.update_process_tensor() + + return self._process_tensor + + +def tti_tempo_compute( + bath: Bath, + start_time: float, + parameters: TempoParameters = None, + progress_type: Optional[Text] = None, + process_tensor_file: Optional[Union[Text, bool]] = None, + overwrite: Optional[bool] = False, + backend_config: Optional[Dict] = None, + name: Optional[Text] = None, + description: Optional[Text] = None) -> BaseProcessTensor: + """ + Shortcut for creating a process tensor by performing a TTI-TEMPO + computation. + + Parameters + ---------- + bath: Bath + The Bath (includes the coupling operator to the system). + start_time: float + The start time. + parameters: TempoParameters + The parameters for the TTI-TEMPO computation. + progress_type: str (default = None) + The progress report type during the computation. Types are: + {``'silent'``, ``'simple'``, ``'bar'``}. If `None` then + the default progress type is used. + name: str (default = None) + An optional name for the tempo object. + description: str (default = None) + An optional description of the tempo object. + """ + ptt = TTITempo(bath, + start_time, + parameters, + process_tensor_file, + overwrite, + backend_config, + name, + description) + ptt.compute(progress_type=progress_type) + return ptt.get_process_tensor() diff --git a/oqupy/util.py b/oqupy/util.py index 07d8c91b..011de820 100644 --- a/oqupy/util.py +++ b/oqupy/util.py @@ -56,6 +56,17 @@ def create_delta( return ret_ndarray +def create_delta_lastindex( + tensor: ndarray) -> ndarray: + """ + Creates a delta at the last index. + """ + ret_shape=tensor.shape+(tensor.shape[-1],) + ret_ndarray=np.zeros(ret_shape,dtype=tensor.dtype) + for a in range(ret_shape[-1]): + ret_ndarray[:,:,a,a]=tensor[:,:,a] + return ret_ndarray + def increase_list_of_index( a: List, shape: List, diff --git a/pylintrc b/pylintrc index 0451ef9c..3607bbf6 100644 --- a/pylintrc +++ b/pylintrc @@ -1,5 +1,4 @@ [MASTER] - # A comma-separated list of package or module names from where C extensions may # be loaded. Extensions are loading into the active Python interpreter and may # run arbitrary code diff --git a/requirements.txt b/requirements.txt index fc7aca66..d3e0bbdd 100644 --- a/requirements.txt +++ b/requirements.txt @@ -4,4 +4,3 @@ tensornetwork==0.4.3 matplotlib>=3.0.0 h5py>=3.1.0 numdifftools>=0.9.41 -