What is the expected behavior?

Once PR #139 is in, add the channel methods to the class, e.g. drive_channel, control_channel, ..., etc. The intended behaviour will need to be determined: e.g. at instantiation should the Solver instance be inspected to see which channels are in the model?

Informations

• Qiskit Ignis version: NA
• Python version:NA
• Operating system:NA

What is the current behavior?

- **Qiskit Ignis version**:? Is this right or it should be Qiskit ODE version or Qiskit terra version

Steps to reproduce the problem

What is the expected behavior?

Suggested solutions

Informations

This is a documentation related issue so no qiskit version info is needed

What is the current behavior?

Tutorials link on top menu is not working, it points to a 404 page

Steps to reproduce the problem

1. Go to https://qiskit.org/documentation/dynamics/
2. Click on Tutorials on top menu

3. Wait until being redirected to the Page not found page

What is the expected behavior?

Instead of having Page not found page, we should get the Tutorials page.

What is the current behavior?

The following raises an error (note that this code depends on syntax from #30):

dispatch.set_default_backend('jax')
static_op = np.array([[1, 0.], [0., -1.]])
operators = [np.array([[0., 1.], [1., 0.]])]

ham_model = HamiltonianModel(static_operator=static_op, 
                             operators=operators,
                             signals=[1.],
                             evaluation_mode='sparse')
ham_model(1.)

This is not a particularly major or urgent bug: sparse operations aren't even supported in JAX anyway (in this package - and in JAX itself there are just experimental implementations of sparse matrices), so nothing is actually being blocked by this.

Nevertheless: I didn't expect this to result in an error. The error is due to the signals being evaluated as an Array if the default backend is JAX. Evaluation of sparse models uses object arrays and np.tensordot, and the Array(..., backend='jax') for the coefficient array causes this to get dispatched to JAX (which does not supported object arrays and therefore raises an error).

Possible solutions

It may be natural to just raise an error if a sparse mode is set while the default backend is JAX. This makes it apparent that there is no JAX sparse support in this package.

the add_sample function needs to be fixed, and tests need to be added to confirm that the padding works.

The code:

discrete = DiscreteSignal(dt=0.5, samples=np.array([]), carrier_freq=3.0)

t_vals = np.array([0.1, .6, 1.1, 1.6])

discrete3.add_samples(1, [4.0, 3.2])

Expected behavior:
discrete.envelope(t_vals) = array([0.0, 4. +0.j, 3.2+0.j, 0. +0.j])

Actual behavior:
discrete.envelope(t_vals) = array([4. +0.j, 3.2+0.j, 0. +0.j, 0. +0.j])

Informations

• Qiskit Dynamics version: main
• Python version: 3.10
• Operating system: Mac OS
• JAX version: 0.3.24

What is the current behavior?

The gradient computation in the how-to-use JAX user guide is returning nan. Weirdly, the gradient computation in the optimization tutorial is still working.

Suggested solutions

This definitely seems to be a Dynamics problem - if I checkout stable/0.3 then it works on the latest version of JAX, with the caveat that another bug (introduced by JAX becoming more strict with type checking) needs to first be fixed to get the code to run.

What is the expected behavior?

Similar to how we are able to use qiskit.quantum_info.Operator types as parameters for models, Operators of base class qiskit.opflow.OperatorBase should also be usable.

Current behaviour

Code Example:

from qiskit_dynamics.solvers import Solver
from qiskit.opflow import Z
drift = Z

hamiltonian_solver = Solver(
    static_hamiltonian=drift
)

Throws DispatchError: 'Dispatch does not support numpy object arrays.'

Solution idea

Add type qiskit.opflow.OperatorBase in function qiskit_dynamics.type_utils.to_numeric_matrix_type

Informations

• Python version:

Python 3.9.15 (main, Oct 11 2022, 21:39:54) 
[Clang 14.0.0 (clang-1400.0.29.102)] on darwin

What is the current behavior?

autosummary is getting confused about the contents of array.Array and failing the build. Previous debugging of this sort of stuff leads me to believe it's only complaining about array because it comes first alphabetically.

Steps to reproduce the problem

Setup:

mkdir tmp && cd tmp
python -m venv .venv
source .venv/bin/activate

git clone [email protected]:Qiskit/qiskit-dynamics.git --depth 1
cd qiskit-dynamics && git rev-parse HEAD
>>> e8a2eb883a9c03385b2cfa98993f051d4d24c829

pip install -r requirements-dev.txt
pip install -e .
pip freeze

alabaster==0.7.12
appnope==0.1.3
astroid==2.9.3
asttokens==2.0.8
attrs==22.1.0
autopage==0.5.1
Babel==2.10.3
backcall==0.2.0
beautifulsoup4==4.11.1
black==22.10.0
bleach==5.0.1
certifi==2022.9.24
charset-normalizer==2.1.1
click==8.1.3
cliff==4.0.0
cmd2==2.4.2
contourpy==1.0.5
cycler==0.11.0
ddt==1.4.4
debugpy==1.6.3
decorator==5.1.1
defusedxml==0.7.1
dill==0.3.5.1
docutils==0.19
dulwich==0.20.46
entrypoints==0.4
executing==1.1.1
extras==1.0.0
fastjsonschema==2.16.2
fixtures==4.0.1
fonttools==4.37.4
future==0.18.2
idna==3.4
imagesize==1.4.1
importlib-metadata==5.0.0
ipykernel==6.16.1
ipython==8.5.0
ipywidgets==8.0.2
isort==5.10.1
jedi==0.18.1
Jinja2==3.1.2
jsonschema==4.16.0
jupyter-sphinx==0.4.0
jupyter_client==7.4.3
jupyter_core==4.11.2
jupyterlab-pygments==0.2.2
jupyterlab-widgets==3.0.3
kiwisolver==1.4.4
latexcodec==2.0.1
lazy-object-proxy==1.7.1
MarkupSafe==2.1.1
matplotlib==3.6.1
matplotlib-inline==0.1.6
mccabe==0.6.1
mistune==2.0.4
mpmath==1.2.1
multiset==3.0.1
mypy-extensions==0.4.3
nbclient==0.7.0
nbconvert==7.2.2
nbformat==5.7.0
nbsphinx==0.8.9
nest-asyncio==1.5.6
numpy==1.23.4
packaging==21.3
pandocfilters==1.5.0
parso==0.8.3
pathspec==0.10.1
pbr==5.11.0
pexpect==4.8.0
pickleshare==0.7.5
Pillow==9.2.0
platformdirs==2.5.2
ply==3.11
prettytable==3.4.1
prompt-toolkit==3.0.31
psutil==5.9.3
ptyprocess==0.7.0
pure-eval==0.2.2
pybtex==0.24.0
pybtex-docutils==1.0.2
Pygments==2.13.0
pylint==2.12.2
pyparsing==3.0.9
pyperclip==1.8.2
pyrsistent==0.18.1
python-dateutil==2.8.2
python-subunit==1.4.0
pytz==2022.5
PyYAML==6.0
pyzmq==24.0.1
-e git+ssh://[email protected]/Qiskit/qiskit-dynamics.git@e8a2eb883a9c03385b2cfa98993f051d4d24c829#egg=qiskit_dynamics
qiskit-sphinx-theme==1.9.0
qiskit-terra==0.22.0
qutip==4.7.0
reno==3.5.0
requests==2.28.1
retworkx==0.12.0
rustworkx==0.12.0
scipy==1.9.3
six==1.16.0
snowballstemmer==2.2.0
soupsieve==2.3.2.post1
Sphinx==5.3.0
sphinx_autodoc_typehints==1.19.4
sphinxcontrib-applehelp==1.0.2
sphinxcontrib-bibtex==2.5.0
sphinxcontrib-devhelp==1.0.2
sphinxcontrib-htmlhelp==2.0.0
sphinxcontrib-jsmath==1.0.1
sphinxcontrib-qthelp==1.0.3
sphinxcontrib-serializinghtml==1.1.5
stack-data==0.5.1
stestr==4.0.1
stevedore==4.1.0
symengine==0.9.2
sympy==1.11.1
testtools==2.5.0
tinycss2==1.2.1
toml==0.10.2
tomli==2.0.1
tornado==6.2
traitlets==5.5.0
typing_extensions==4.4.0
urllib3==1.26.12
voluptuous==0.13.1
wcwidth==0.2.5
webencodings==0.5.1
widgetsnbextension==4.0.3
wrapt==1.13.3
zipp==3.9.0

Build attempt:

python -m sphinx -b html -W docs/ docs/_build/html

Returns:

Running Sphinx v5.3.0
making output directory... done
checking bibtex cache... out of date
parsing bibtex file /Users/ihincks/ttmp/qiskit-dynamics/docs/refs.bib... parsed 12 entries
[autosummary] generating autosummary for: apidocs/array.rst, apidocs/dispatch.rst, apidocs/index.rst, apidocs/models.rst, apidocs/perturbation.rst, apidocs/pulse.rst, apidocs/signals.rst, apidocs/solvers.rst, discussions/dyson_magnus.rst, discussions/index.rst, ..., release_notes.rst, tutorials/Lindblad_dynamics_simulation.rst, tutorials/Rabi_oscillations.rst, tutorials/index.rst, tutorials/optimizing_pulse_sequence.rst, tutorials/qiskit_pulse.rst, userguide/how_to_configure_simulations.rst, userguide/how_to_use_jax.rst, userguide/index.rst, userguide/perturbative_solvers.rst
[autosummary] generating autosummary for: /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.array.Array.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.array.wrap.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.dispatch.asarray.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.dispatch.requires_backend.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.models.GeneratorModel.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.models.HamiltonianModel.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.models.LindbladModel.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.models.RotatingFrame.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.models.rotating_wave_approximation.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.perturbation.ArrayPolynomial.rst, ..., /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.signals.DiscreteSignalSum.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.signals.IQMixer.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.signals.Signal.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.signals.SignalList.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.signals.SignalSum.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.solvers.DysonSolver.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.solvers.MagnusSolver.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.solvers.Solver.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.solvers.solve_lmde.rst, /Users/ihincks/ttmp/qiskit-dynamics/docs/stubs/qiskit_dynamics.solvers.solve_ode.rst

Warning, treated as error:
[autosummary] failed to import qiskit_dynamics.array.Array.available_backends.
Possible hints:
* ImportError: 
* AttributeError: module 'qiskit_dynamics.array.Array' has no attribute 'available_backends'
* ModuleNotFoundError: No module named 'qiskit_dynamics.array.Array.available_backends'; 'qiskit_dynamics.array.Array' is not a package
*

As usual, the -P option on sphinx is useless.

Unsurprisingly, same behaviour with tox as the environment is basically identical:

pip intsall tox
(.venv) ihincks@Ians-MacBook-Pro qiskit-dynamics % which tox
/Users/ihincks/ttmp/.venv/bin/tox
(.venv) ihincks@Ians-MacBook-Pro qiskit-dynamics % tox -edocs
docs create: /Users/ihincks/ttmp/qiskit-dynamics/.tox/docs
docs installdeps: -r/Users/ihincks/ttmp/qiskit-dynamics/requirements-dev.txt, jax, jaxlib, diffrax
docs develop-inst: /Users/ihincks/ttmp/qiskit-dynamics
docs installed: absl-py==1.3.0,alabaster==0.7.12,appnope==0.1.3,astroid==2.9.3,asttokens==2.0.8,attrs==22.1.0,autopage==0.5.1,Babel==2.10.3,backcall==0.2.0,beautifulsoup4==4.11.1,black==22.10.0,bleach==5.0.1,certifi==2022.9.24,charset-normalizer==2.1.1,click==8.1.3,cliff==4.0.0,cmd2==2.4.2,contourpy==1.0.5,cycler==0.11.0,ddt==1.4.4,debugpy==1.6.3,decorator==5.1.1,defusedxml==0.7.1,diffrax==0.2.1,dill==0.3.5.1,docutils==0.19,dulwich==0.20.46,entrypoints==0.4,equinox==0.8.0,etils==0.8.0,executing==1.1.1,extras==1.0.0,fastjsonschema==2.16.2,fixtures==4.0.1,fonttools==4.37.4,future==0.18.2,idna==3.4,imagesize==1.4.1,importlib-metadata==5.0.0,importlib-resources==5.10.0,ipykernel==6.16.1,ipython==8.5.0,ipywidgets==8.0.2,isort==5.10.1,jax==0.3.23,jaxlib==0.3.22,jaxtyping==0.2.7,jedi==0.18.1,Jinja2==3.1.2,jsonschema==4.16.0,jupyter-sphinx==0.4.0,jupyter_client==7.4.3,jupyter_core==4.11.2,jupyterlab-pygments==0.2.2,jupyterlab-widgets==3.0.3,kiwisolver==1.4.4,latexcodec==2.0.1,lazy-object-proxy==1.7.1,MarkupSafe==2.1.1,matplotlib==3.6.1,matplotlib-inline==0.1.6,mccabe==0.6.1,mistune==2.0.4,mpmath==1.2.1,multiset==3.0.1,mypy-extensions==0.4.3,nbclient==0.7.0,nbconvert==7.2.2,nbformat==5.7.0,nbsphinx==0.8.9,nest-asyncio==1.5.6,numpy==1.23.4,opt-einsum==3.3.0,packaging==21.3,pandocfilters==1.5.0,parso==0.8.3,pathspec==0.10.1,pbr==5.11.0,pexpect==4.8.0,pickleshare==0.7.5,Pillow==9.2.0,platformdirs==2.5.2,ply==3.11,prettytable==3.4.1,prompt-toolkit==3.0.31,psutil==5.9.3,ptyprocess==0.7.0,pure-eval==0.2.2,pybtex==0.24.0,pybtex-docutils==1.0.2,Pygments==2.13.0,pylint==2.12.2,pyparsing==3.0.9,pyperclip==1.8.2,pyrsistent==0.18.1,python-dateutil==2.8.2,python-subunit==1.4.0,pytz==2022.5,PyYAML==6.0,pyzmq==24.0.1,-e git+ssh://[email protected]/Qiskit/qiskit-dynamics.git@e8a2eb883a9c03385b2cfa98993f051d4d24c829#egg=qiskit_dynamics,qiskit-sphinx-theme==1.9.0,qiskit-terra==0.22.0,qutip==4.7.0,reno==3.5.0,requests==2.28.1,retworkx==0.12.0,rustworkx==0.12.0,scipy==1.9.3,six==1.16.0,snowballstemmer==2.2.0,soupsieve==2.3.2.post1,Sphinx==5.3.0,sphinx_autodoc_typehints==1.19.4,sphinxcontrib-applehelp==1.0.2,sphinxcontrib-bibtex==2.5.0,sphinxcontrib-devhelp==1.0.2,sphinxcontrib-htmlhelp==2.0.0,sphinxcontrib-jsmath==1.0.1,sphinxcontrib-qthelp==1.0.3,sphinxcontrib-serializinghtml==1.1.5,stack-data==0.5.1,stestr==4.0.1,stevedore==4.1.0,symengine==0.9.2,sympy==1.11.1,testtools==2.5.0,tinycss2==1.2.1,toml==0.10.2,tomli==2.0.1,tornado==6.2,traitlets==5.5.0,typeguard==2.13.3,typing_extensions==4.4.0,urllib3==1.26.12,voluptuous==0.13.1,wcwidth==0.2.5,webencodings==0.5.1,widgetsnbextension==4.0.3,wrapt==1.13.3,zipp==3.9.0
docs run-test-pre: PYTHONHASHSEED='1448775009'
docs run-test: commands[0] | sphinx-build -b html -W docs/ docs/_build/html
Running Sphinx v5.3.0
checking bibtex cache... out of date
parsing bibtex file /Users/ihincks/ttmp/qiskit-dynamics/docs/refs.bib... parsed 12 entries
[autosummary] generating autosummary for: apidocs/array.rst, apidocs/dispatch.rst, apidocs/index.rst, apidocs/models.rst, apidocs/perturbation.rst, apidocs/pulse.rst, apidocs/signals.rst, apidocs/solvers.rst, discussions/dyson_magnus.rst, discussions/index.rst, ..., stubs/qiskit_dynamics.solvers.solve_ode.rst, tutorials/Lindblad_dynamics_simulation.rst, tutorials/Rabi_oscillations.rst, tutorials/index.rst, tutorials/optimizing_pulse_sequence.rst, tutorials/qiskit_pulse.rst, userguide/how_to_configure_simulations.rst, userguide/how_to_use_jax.rst, userguide/index.rst, userguide/perturbative_solvers.rst

Warning, treated as error:
[autosummary] failed to import qiskit_dynamics.array.Array.available_backends.
Possible hints:
* ModuleNotFoundError: No module named 'qiskit_dynamics.array.Array.available_backends'; 'qiskit_dynamics.array.Array' is not a package
* ImportError: 
* AttributeError: module 'qiskit_dynamics.array.Array' has no attribute 'available_backends'
ERROR: InvocationError for command /Users/ihincks/ttmp/qiskit-dynamics/.tox/docs/bin/sphinx-build -b html -W docs/ docs/_build/html (exited with code 2)
_________________________________________________________ summary _________________________________________________________
ERROR:   docs: commands failed

Gives

What is the expected behavior?

Docs should build.

Suggested solutions

What is the expected behavior?

Connected to #53.

Update the converter to include optional argument channels, specifying a list of channels as strings e.g. ['D0', 'D2', 'U0'], which modifies the behaviour of the function to:

• Filter out channels not included in the list.
• Return signals in a list given by the ordering in channels.

Currently the interface for setting signals of HamiltonianModel and LindbladModel uses the setter for the .signals property. This is a bit awkward in the case of the LindbladModel, which has two collections of signals for Hamiltonian and dissipator terms. A result is that, the signals property for Lindblad is now a Tuple of these two terms, and setting it requires setting both, though this is a bit awkward given that a huge use case for the Lindblad model will be when the signals are constantly 1.. It may be good in general to change the signal setter for GeneratorModel/HamiltonianModel to a method set_signals, and Lindblad model should maybe have two separate properties/setters for hamiltonian_signals and dissipator_signals.

What is the expected behavior?

Once PR #139 is merged, add functionality for setting meas_level==1 to the PulseSimulator class options. Should also add meas_return option with a default (I believe this impacts how outcomes are returned when meas_level==1).

For now this can enable IQ data spoofing, and later on we can add the ability to specify custom IQ-data generators.

What is the expected behavior?

In discussion #95 it is pointed out that the unwrapping/wrapping of input states by Solver.solve doesn't preserve metadata of the state. E.g. an input given by a DensityMatrix with subsystem dimensions defined will result in an output without subsystem dimensions.

What is the expected behavior?

Taking/modifying code from Aer, write a function with signature:

• Inputs:
  • Hamiltonian string representation
  • Subsystem list (list of subsystems in the model to keep when parsing)
• Outputs:
  • Hamiltonian static_operator and operators (i.e. operator data required to instantiate a HamiltonianModel or Solver)
  • List of channel names (as strings, e.g. ['D0', 'D2', 'U0'] corresponding to the operators in operators.

The schedule -> signal converter currently ensures that the output signals satisfy Equation (1) of this paper, with the frequency interpreted as the current channel frequency value given by SetFrequency/ShiftFrequency instructions, and the phase interpreted as being the value given purely by SetPhase/ShiftPhase instructions.

As pointed out by @nkanazawa1989 in this notebook however, frequency changes result in phase accumulation on real IBM backends. Hence, the phase in Equation (1) should actually be a function of the SetPhase/ShiftPhase instructions, as well as the SetFrequency/ShiftFrequency instructions (and the times at which they occur).

We should update the behaviour of the pulse -> signal converter to be consistent with this.

Solver is emerging as the main higher level interface we are promoting, but the module API docs describing it in solvers.__init__ are pretty sparse. They don't describe how the internal models are being constructed, or how they're being transformed by rotating frames or the RWA.

All of this information is currently contained in the models.__init__ file. The original idea was that a user could go there to read about this, but with the Solver API becoming the dominant mode of interaction with the package it makes sense to either reproduce this in solvers.__init__, move it there, or something else.

Several small changes needed:

• Change diffrax solver handling of t_eval so that it is not used in the final call to odeResult, causing an issue if t_eval is None
• Add handling for the diffrax solver in solve_lmde
• Add tests for solve_lmde's delegation to solve_ode, basic cases for each solver to confirm the pipeline works
• Add diffrax solver to the solve_ode docs

Currently, the hamiltonian operator validation step will raise an error if attempting to JAX-transform a function in which a hamiltonian/lindblad model or a Solver is constructed. Setting validate=False bypasses this, and in some cases leads to being able to JAX-transform a function in which one of these objects is constructed. However, some cases that should work don't, e.g. passing in operators=[Operator.from_label('X')] to the constructor leads to an error, which seems to stem from a tracing issue with to_array.

An example of a case that works is in test.dynamics.solvers.test_solver_classes.TestSolverJax.test_transform_through_construction_when_validate_false (currently in PR #30 but should be merged soon). We should add more tests of this form with different input types, or identify what particular part of construction is causing an error and add jit tests for that (e.g. if to_array is the problem, add tests for using to_array within a function and try transforming it).

The following code results in an error:

import numpy as np

from qiskit_dynamics import Solver
from qiskit_dynamics import dispatch

solver = Solver(hamiltonian_operators=[np.array([[0., 1.], [1., 0.]])], 
                hamiltonian_signals=[1.])

dispatch.set_default_backend('jax')

solver.solve(t_span=[0., 1.], y0=np.array([0., 1.]), method='jax_odeint')

Moving dispatch.set_default_backend('jax') before the Solver construction results in the above working.

The error comes from internal choices being made in Solver based on what the default backend is at the time of instantiation. E.g.:

• The way in which signals are evaluated changes depending on whether the backend is JAX or numpy.
• The backend for the operators stored in the internal models is also set at instantiation, and changing the default backend after instantiation doesn't change that.
• There may be others.

I'm labelling this as a bug as it results in an error, but I mainly view this as an "enhancement" in the sense that the intended usage is to set the default backend at the beginning, and this is outside of that usage.

What is the expected behavior?

When taking the matrix exponential of an Hermitean matrix $M$ (for example to compute the time evolution operator from the Hamiltonian of a closed system) it is sometimes faster to first diagonalize the matrix $M = V D V^{dagger}$ and then exponentiate the diagonal matrix $D$ (which is given by exponentiating the diagonal elements) instead of using the built-in method in scipy expm directly on the matrix $M$.

Since fixed-step solvers are meant to work for arbitrary generators (even ones that cannot be diagonalized) but this method is faster for diagonalizable generators (in particular anti-hermitian), then it would be worth adding this implementation as an option for the user and to choose this option automatically in some higher level interfaces (e.g. when using a Solver instance that only has Hamiltonian information).

Possible implementation details:
Add the option to the core solver method, and add a check in the solve_lmde interface that if what's being simulated is a HamiltonianModel object, and the dim of the initial state is below some heuristically-set level, then it automatically uses the diagonalization method.

What is the expected behavior?

I suggest the builder classes (in the context of Hamiltonian/Lindblad models) to support the following workflow:

1. Define Hamiltonian/Lindblad operators.
2. Do (potentially costly) precomputations (build operators), targeting different solvers and/or varying build parameters while using a single definition of the model operators from step (a).
3. Specify an initial state, simulation times, solver parameters, and solve, possibly multiple times (varying input), using one set of prebuilt operators from step (b).
4. Potentially serialize a solution from step (c) to files.
5. Potentially load a serialized solution from files.
6. Specify and calculate observables using a solution.
7. Extract observables from multiple (different) simulations, do simple manipulations on the data, and plot.

Builders should also balance code reuse of similar tasks between different solvers. Specifically, try to encompass the important similarities when solving Hamiltonian/Lindblad models (despite differences in details between solvers).

When benchmarking sparse computations on jax vs qutip, a slowdown was observed when the Signal class was used in the qutip simulation. Replacing the Signal usage with a manually defined numpy pulse function gave a 6x speedup for a small 2-transmon simulation with 3-d qubits. In addition, vectorization also caused a slowdown, with the np.vectorization wrapper causing a 2x slowdown in the same process. This testing was performed with qutip, using the qt.mesolve function.

Currently DiscreteSignal.envelope is setup to clip the time into the time-range on which the samples are defined. The result is that for all time before start time, the envelope evaluates to the first sample, and for all times after the end time, it evaluates to the last sample.

E.g. For a discrete signal with start_time=0, samples=[1., 2., 3.] and dt=1., calling discrete_signal.envelope(4.) will return 3., and discrete_signal.envelope(-1.) will return 1..

It will make more sense, and be more useful for coming applications, if the envelope returns 0. when t is outside the sample window.

Informations

• Qiskit Dynamics version: stable/0.3
• Python version: 3.10
• Operating system: Mac OS

What is the current behavior?

The how-to-use JAX userguide raises an error related to jnp.ones_like not accepting Array as an input. This was resolved in PR #130, but we should do a minor or patch release to 0.3 to fix this, either putting an upper bound on the JAX version, or by including the PRs adding fixes for JAX updates.

What is the expected behavior?

Add a function (or enhance an existing one) to RotatingFrame that is qiskit.quantum_info "type aware", to be used at a higher level so that a user doesn't need to change syntax based on type.

Currently RotatingFrame has a bunch of frame-transformation helper functions that explicitly perform certain mathematical operations, e.g.:

rotating_frame.state_into_frame(t, y)

performs the mapping t, y -> exp(-tF) @ y and

rotating_frame.operator_into_frame(t, y)

performs the mapping t, y -> exp(-tF) @ y @ exp(tF) (might have minus signs mixed up).

Currently, if a user has a Statevector instance, they need to call state_into_frame, but if they have a DensityMatrix, they need to call operator_into_frame to perform the correct transformation. Furthermore, all current functions will strip any wrapping, and return a raw array. Note as well that the terminology of these functions is meant more to differentiate between the mathematical transformations, rather than to align with standard quantum usage (e.g. a DensityMatrix is a "state").

It'd be great to have something higher level, like rotating_frame.into_frame(t, y) (and similarly rotating_frame.out_of_frame(t, y), which utilizes the type of y to determine the correct mathematical transformation, and returns the results in the original wrapped class.

As pointed out by @chriseclectic , the function InstructionToSignals.get_signals has been updated so that the carriers argument is supplied as a dictionary mapping channel names to frequencies, rather than a list of frequencies (whose channel associations are automatically determined by the iteration ordering of schedule.channels). This is a breaking change (and was discovered by running main code while playing around with the existing pulse to signals tutorial).

We should handle this change with a deprecation warning:

• Re-add the ability to pass carriers as a list by creating the proper dict format using the associations the old code would have generated.
• Include a deprecation warning, that this style of formatting is deprecated as of 0.3.0, and will be fully removed by 0.4.0.
• Update the release note file for the pulse-converter-filtering to include a note that this deprecation is occurring.

What is the expected behavior?

Compute probabilities function:

Inputs:

• np.ndarray (1d for state vector or 2d for density matrix) or one of the quantum info QuantumState types.
• 'dressed states' (eigenvectors for drift hamiltonian that defines the measurement projectors) (output of function in #74 )
• 'dressed state labels' (the computational state strings corresponding to the dressed states) (output of function in #74 )

Outputs:

• Vector of probabilities (ordered appropriately according to 'dressed state labels')

Sample counts:

Inputs:

• Probability vector with corresponding outcome labels
• n_shots - number of samples to draw

Output:

• Dict of counts of samples, with keys being outcome labels and values being the counts

What is the expected behavior?

Introduce higher level solver classes for solving quantum systems that assemble/utilize lower level components. This is the beginning of an interface layer that will provide some uniformity across different execution paths (e.g. what is in the package currently v.s. MPO v.s. perturbative v.s. unknown) - while these underlying methods are different and therefore the higher level interfaces here may need to be different, they can at least be somewhat consolidated.

In conjunction with this, a lot of the interface flair of solve_lmde can be removed and inserted into the higher level solver interface. This function needs to be simplified in any case, as the logic for setting up the simulation based on optional arguments is unnecessarily convoluted and therefore confusing to maintain.

Initial plan

We can create a stateful solver class that can be updated and stepped forward. It can internally set up HamiltonianModel or LindbladModel as required, and handle a lot of the high level interface considerations originally present in solve_lmde (e.g. automatic wrapping/unwrapping of terra state classes, etc...).

One additional consideration is how to handle JAX execution, as JAX compilation causes problems if part of the function being compiled involves setting properties of objects. Some options here:

• We can make a JaxSolver which, after initialization, contains only pure methods. The benefit of this is that it'll make it very clear what "can be compiled" v.s. what can't, though at the price of defining a separate interface.
• We can advise a usage mode which is similar to what currently needs to be done with models: if a function to be compiled involves setting a property of the Solver (either, e.g., directly setting signals, or simulating and storing the latest y value), that function can be "purified" by first making a copy of the Solver. This will allow us to use a single class/interface, at the cost of adding a principle for "standard JAX usage" (though this principle essentially already exists with the models).

Things that need be changed, as provided by @jakelishman

• tox.ini: remove the py36 environment
• setup.py: remove the “Python 3.6” PyPI classifier and update python_requires
  • Looks like we need to actually add "Python 3.9" to this as well.
• pyproject.toml: update black’s configuration
• .github/workflows/*.yml: update any python-version matrix specifiers
• CONTRIBUTING.md: remove references to py34 and py35

Passing an empty DiscreteSignal into to_SignalSum causes an array with dimensions (1,0) to be transpoed to (0,1), which breaks the DiscreteSignal init function. The transposition should not be applied in the empty sample case.

Current behavior (raises error):

from qiskit_dynamics import DiscreteSignal
from qiskit_dynamics.signals.signals import to_SignalSum

to_SignalSum(DiscreteSignal(samples=[], dt=1.))

To fix, add an if statement to to_SignalSum which only transposes the input samples if their size is not (1,0)

A new user recently had a problem with running tutorial code (both of the first two tutorials numpy/scipy based simulations), getting the error

TypeError: ufunc 'nextafter' not supported for the input types, and the inputs could not be safely coerced to any supported types according to the casting rule ''safe''

Their numpy version was 1.19.2 (which is well above the minimum version listed for dynamics), but upgrading resolved the issue. We should find the true minimum and update the listing to avoid this problem.

Figure out a minimal JAX version that seems to work and add CI tests for this.

During development we have been writing documentation and example jupyter notebooks, and for the first release we will need to polish the docs and move the example notebooks to tutorials/docs as per general qiskit style.

Based on discussion, some rough organizing principles that might be good are (these are not hard rules, but can guide what goes where):

• Tutorials should be holistic demonstrations of a use-case of the package, and should not be focused on outlining the interface/functionality of a specific class or function (which are the domain of docs).
• Module-level docs demonstrate broadly what the sub-module has to offer, with a focus on usability/functionality. As they are less concerned with source readability than class/method/function docs (see below), they can make more liberal usage of things like jupyter-execute for displaying plots, etc.
• Class/method/function docs, aside from obviously being very specific to the given object they are describing, should aim to maintain legibility of help calls on the object.

The existing example notebooks are a mix of "module-level documentation" and "tutorials", and are probably even leaning more towards the former.

To do:

• Update readme.
• Add docs/tutorial folder and CI tests for tutorials. (PR #16)
• Make a specific plan as final major code contributions are finishing up.
  • Final decision on organizing principles.
  • Decide on what we generally want out of tutorials/docs.
• Implement plan!
  • Migrate existing examples notebooks into proper doc/tutorial locations.
  • Add new tutorials/docs as necessary.

What is the expected behavior?

There are a few type-handling issues that we need to sort out. Currently, the model classes support operation in both dense and sparse modes, however, regardless of the chosen mode of operation, at instantiation, everything is passed through the to_array function, which converts everything to a dense array. If a sparse mode is selected, it will then convert these arrays to sparse. We should set things up so that if a user specifies things in a sparse format then it never gets converted to a dense format.

To do:

• Add a to_csr analog of to_array, which does the same job but converts everything to the required format for the sparse types (see e.g. the constructor to SparseOperatorCollection for the required format).
• Add handling of qutip operator input types to both to_array and to_csr. This should be done without actually importing qutip, so need to decide how to detect this (op.__name__ ?) .
• Make sure the use-case of specifying things operators in a sparse format, and requesting models to operate in sparse format, never converts to dense.

Informations

This is a documentation related issue so no qiskit version info is needed

What is the current behavior?

Getting started link on top menu is not working, it points to a 404 page

Steps to reproduce the problem

1. Go to https://qiskit.org/documentation/dynamics/
2. Click on Getting started on top menu

3. Wait until being redirected to the Page not found page

What is the expected behavior?

Instead of having Page not found page, we should have a Getting started page. If we don't have a Getting started page, we should remove the option from the top menu

PR #51 introduces JAX sparse evaluation for models. Due to some limited functionality for the jax BCOO sparse array type, some work arounds were required to implement a few things. Subsequent releases of JAX are expected to eliminate the need for these workarounds, and this issue is a reminder of what these are.

• A function jsparse_linear_combo is defined at the beginning of operator_collections.py that takes a linear combination of sparse arrays (specified as a 3d BCOO array), with the coefficients given in a dense 1d array. This cannot be achieved using a sparisfied version of jnp.tensordot as the design convention jax is using is that such operations will always output dense arrays if at least one input is dense. Hence, jsparse_linear_combo multiplies the coefficients against the sparse array directly via broadcasting. However, sparse-dense element-wise multiplication, at the time of writing, is limited to arrays of the same shape, and therefore it is necessary to explicitly blow up the coefficient array to a dense array with the same shape as the sparse array (which is huge). I'm not sure if this is done via views so it's okay, but this should be changed when possible regardless.
• Setting up of operators in the jax sparse Lindblad operator collection is done using dense arrays, as sparse-sparse matrix multiplication is not yet implemented. This is relatively minor but it would be nice to do it with sparse when possible.
• "Vectorized" products A @ B, where A and B are 3d arrays, with one being sparse and the other dense, are not reverse-mode differentiable. This results in LindbladModel in jax-sparse not being reverse-mode differentiable. It is however, forward mode differentiable. At some point this will change, and we will need to remove the caveat in LindbladModel.evaluation_mode that sparse mode with jax is not reverse-mode differentiable.

When simulating a list of pulse schedules using JAX, the Solver.solve method internally takes advantage of JAX compilation to speed up simulation. To get around the limitation of compilation that require function inputs to be a fixed shape, the strategy is as follows:

• Given a list of schedules, determine the max duration.
• Build a function for simulating a single schedule with length of the max duration.
  • The inputs of this function include a large array containing all samples for all channels in the schedule, along with the duration of the specific schedule in units of time T.
• Tell JAX to compile the function.
• For each schedule, extract all of the samples into a single array, padding them to be the length of the max duration schedule (so they are all the same), and call the simulation function with this array and the integration interval [0, T].

As all samples are always padded out to be the length of the max duration schedule, JAX will only compile the function on the first schedule, then re-use the compiled function for all remaining schedules.

This strategy works well to speed up a single call to Solver.solve, but doesn't do anything to help repeated calls: the compilation cost is incurred with every call to Solver.solve. For small system sizes, where compilation can literally take orders of magnitude longer than simulation itself (e.g. 1 second vs 10 ms), this cost greatly impacts the speed of JAX-based pulse simulation in Dynamics.

As discussed with @ihincks , alternative strategies could fix this issue. For example:

• Compile the simulation function to a fixed size.
• For each encountered schedule: if it is a shorter duration, pad it and simulate with the compiled function, and if it is a longer duration, recompile the function for a longer fixed size (chosen heuristically).
• Repeat, possibly keeping track of each duration that was used for compilation, and always padding each schedule with the shortest number of samples to achieve an already-compiled-over shape.

Some version of this would be a great update.

The old general_demo.ipynb did a good job of laying out the pieces/core capabilities of the package. However, it ended up being deleted as:

• It worked directly with models and solve_lmde, and this is now out of date since Solver was introduced.
• It didn't really fit into any documentation category: it was a mix of tutorial/userguide/discussion.

I think though that we really need an updated version of this somewhere in the documentation. Probably either a User guide entry or a Discussion piece makes sense. I'm leaning towards "discussion" piece as I'm still conceptualizing "user guide" to be dedicated to more advanced/technical information about how to use the package in question/answer format, whereas this is more of a general description of how the package is organized.

I just spent several hours being confused by a bug that could have been avoided if an incompatible shape error was raised when constructing instances of the DysonSolver and MagnusSolver classes.

Specifically, I had a model with 4 operators, but when I built the solver, I accidentally only specified a list with two ints for chebyshev_orders, rather than 4 (to match the length of operators). The instance was still built, and successfully ran simulations, but the signals beyond the first two were always ignored (i.e. the simulation ran as if those two signals were 0. at all times).

It's kind of interesting that it still ran - but in any case, DysonSolver, MagnusSolver, and the underlying ExpansionModel definitely implicitly assume that operators, carrier_freqs, and chebyshev_orders all have the same length. Validation checks should be added to ExpansionModel that verify this.

Summary

Make the fixed step solvers compileable/differentiable with respect to the t_span and t_eval arguments.

This isn't an extremely urgent issue, however it would still be very nice to round out the features of these solvers.

Details

Issue #122 outlines a bug that is ultimately due to the fact that the JAX solvers in dynamics cannot be compiled if t_eval is not None. The fix PR, #125, resolves this issue by updating jax_odeint and the diffrax solver wrapper so that they can be compiled with respect to t_eval.

As described in #122 however, updating the fixed step JAX solvers built in dynamics to be compilable with respect to both t_span and t_eval is non-trivial due to their looping structure being dependent on the values of t_span and t_eval. As a result, the fix #125 is only partial: in the case of JAX fixed step solvers, the problem is simply avoided rather than being fundamentally fixed.

To make the fixed step solvers fully compilable/differentiable with respect to the t_span and t_eval arguments, the functions fixed_step_solver_template_jax and fixed_step_lmde_solver_parallel_template_jax need to be updated to use more advanced JAX control flow. Preserving differentiability with respect to other parameters may also require defining custom differentiation rules (vjp and jvp rules).

A group of small API changes:

• Argument naming conventions
  • Change drift kwarg to static_hamiltonian or static_operator where applicable
• Somewhat related to the above:
  • Add static dissipator for lindblad?
  • All model classes should be able to be instantiated with/without various operator types. E.g. currently operators is required for HamiltonianModel, but we should be able to instantiate with just a drift. Similarly, a LindbladModel should be able to be instantiated with just dissipators (no Hamiltonian part).
    • It should be "simple" to change this but may be somewhat non-trivial as the current structure implicitly makes these assumptions in a bunch of places.
• Add validation to Lindblad signal setter
  • None currently exists.
• Add validation to Hamiltonian operators in Lindblad
  • None currently exists.
• For some reason Lindblad constructor needs Hamiltonian signals.
  • This should not be necessary.

Depending on the scope of these they could be handled in a single PR, but if any require substantial changes they should be done in subsequent PRs.

A note should be added to the userguide (using the sphinx .. note:: environment so it stands out) highlighting that JAX defaults to operate in 32 bit mode, and that it needs to explicitly be set to 64 bit mode to get high accuracy solutions.

The userguide currently shows how to set to 64 bit mode, but doesn't emphasize why this is necessary.

With #164 , we're working towards making the InstructionToSignals class behave in a way more reflective of how schedules are implemented on real backends. As this behaviour is non-obvious and non-trivial, after #164, we should update the documentation for this class to be much more explicit about the rules its following to generate the samples used in simulation from schedules.

What is the expected behavior?

When a parallel solver is used, the solver attempts to parallelize over every single time step that is being solved. As a result, these parallel solvers use massive amounts of memory. In order to enable these solvers to be used at higher dimensions, or over longer time periods, some functionality should be added which limits how many steps are parallelized at once.

A couple different ways to do this:

• Have the user manually input how many steps to parallelize over -- fast solution but less user-friendly as it requires a lot of minute control, may be a good initial solve
• Have the user input how much memory they have available and then have a heuristic to choose parallelization based on that -- users must note that other parts of dynamics may use additional memory, and it may be hard to calculate exactly how much memory the parallel solver uses
• Fully automate the process based on some internal Jax memory limits

What is the expected behavior?

Inputs:

• Static hamiltonian - a hermitian operator that is assumed to be nearly diagonal
• Subsystem dimensions - a list of the dimensions of the tensor-factors

Outputs:

• The 'dressed states': eigenvectors of the static hamiltonian labelled according to which computational state they correspond to.
  • E.g. this could be a dictionary whose keys are the computational state labels (strings of 0, 1, 2, ...) with values being the corresponding eigenstate
  • Alternatively, same as the above but given as two lists.
• The 'dressed frequencies': For each subsystem, the difference between the eigenvalues of the dressed first excited state and the ground state, divided by 2 pi. This could just be returned as a 1d array.

Some code for performing these steps already exists in Aer in the HamiltonianModel._compute_drift_data() method, so we can draw on this if it makes sense to do so.

Informations

• Qiskit Dynamics version: 0.3.0
• Python version: 3.10
• Operating system: Mac OS X

What is the current behavior?

The automatic jit compilation routine for simulating pulse schedules with JAX raises errors under certain conditions.

Steps to reproduce the problem

Specifying if using method='jax_odeint', specifying the t_eval kwarg will cause the error to be raised. Also, if using method='jax_expm', the error will be raised regardless of whether t_eval is included or not.

Examples here can probably be reduced further in terms of lines of code, but here is a minimal example:

import numpy as np
from qiskit import pulse

from qiskit_dynamics.array import Array
from qiskit_dynamics import Solver

import jax
jax.config.update("jax_enable_x64", True)

# tell JAX we are using CPU
jax.config.update('jax_platform_name', 'cpu')

from qiskit_dynamics.array import Array
Array.set_default_backend('jax')

from qiskit_dynamics import Solver

Z = np.array([[1., 0.], [0., -1.]])
X = np.array([[0., 1.], [1., 0.]])

v = 5.
anharm = -0.33
r = 0.02
dt = 0.222222222222222 # dt for pulse schedules

# static part
static_hamiltonian = 2 * np.pi * v * Z
# drive term
drive_hamiltonian = 2 * np.pi * r * X

solver = Solver(
    static_hamiltonian=static_hamiltonian,
    hamiltonian_operators=[drive_hamiltonian],
    rotating_frame=static_hamiltonian,
    hamiltonian_channels=['d0'],
    channel_carrier_freqs={'d0': v},
    dt=dt
)

with pulse.build(name='x') as schedule:
    pulse.set_frequency(v, channel=pulse.DriveChannel(0))
    pulse.play(pulse.Drag(300, 1., 100, 0.), pulse.DriveChannel(0))

This works:

solver.solve(
    signals=schedule,
    y0=np.eye(2, dtype=complex),
    t_span=[0, 300 * dt],
    method='jax_odeint',
    atol=1e-10,
    rtol=1e-10
)

However this does not:

solver.solve(
    signals=schedule,
    y0=np.eye(2, dtype=complex),
    t_span=[0, 300 * dt],
    t_eval=[0, 150 * dt, 300 * dt],
    method='jax_odeint',
    atol=1e-10,
    rtol=1e-10
)

and similarly this does not:

solver.solve(
    signals=schedule,
    y0=np.eye(2, dtype=complex),
    t_span=[0, 300 * dt],
    t_eval=[0, 150 * dt, 300 * dt],
    method='jax_expm',
    max_dt=dt
)

What is the expected behavior?

The above calls should work and produce correct results.

Suggested solutions

Both errors are being raised in different places, but the underlying reason is the same: the function being internally compiled (located here) takes t_span as an argument, but in both cases (if t_eval is specified, or if a fixed step solver is being used), t_span cannot be compiled over.

Some imperfect options:

• Set t_span as a static argument in the internally compiled function. This is an extremely easy option, but will limit the usefulness of the compilation to cases where all t_span values are the same. It will still apply to a lot of cases, but I think this is an option of last resort.
• Figure out how to make both of the failure cases compilable. I think the first one is possible - the only reason it wasn't written to be compilable originally was due to difficulty, but I think in principle it can be done. This is an imperfect option however as I seem to remember there being fundamental issues with making the second one compilable (will need to see).

What is the current behavior?

Calling solve_lmde with a HamiltonianModel and DensityMatrix initial condition results in solving the DE p'(t) = -iH(t)p(t) (where p(t) is the density matrix). In a sense this is intended behaviour: if the user passes an initial condition for which the generator-state matrix product G(t)y0 makes sense, it will solve y'(t) = G(t)y(t). However, in this case, since the initial condition is wrapped in DensityMatrix, this probably isn't what the user means to solve.

Suggested solutions

As suggested by @haggaila , we could raise an error here, instructing the user to put the HamiltonianModel into a LindbladModel, to be explicit that they want to simulate the Lindblad equation.

What is the expected behavior?

Currently, if a Solver is instantiated with Lindblad data, calling Solver.solve with y0 being a SuperOp will only work if evaluation_mode is set to a vectorized option, and otherwise it will raise an error explaining this.

The internal behaviour of Solver.solve should be updated so that the simulation is still carried forward even if evaluation_mode is not vectorized by:

• Converting the SuperOp to a 3d array representing a list of input matrices to the un-vectorized Lindblad equation.
• Passing this as the initial state to the call to solve_lmde
• Converting the output back to a SuperOp by inverting the first transformation.

Notes

The second bullet point above should work due to the LindbladModel.evaluate_rhs being "vectorized" in the array-programming sense: evaluating it on a 3d array is functionally equivalent to looping over the first index and applying the RHS function to each entry.

• Part of this PR will involve more carefully testing/verifying the components that go into this. There don't seem to be any tests for this in test_lindblad_model.py, but the test test_multiple_density_matrix_evaluation in test.dynamics.models.test_operator_collections.TestDenseLindbladCollection tests this at the operator collection level. The LindbladModel array-vectorized rhs operation should therefore work assuming that the RotatingFrame operator frame transformation functions are also vectorized (which I think they are).

Additionally, once working, it would be good to test out some examples comparing vectorized evaluation_mode v.s. non-vectorized, to get a sense of the performance difference (though this isn't technically required for this).

PR #51 introduces JAX sparse evaluation mode. As suggested by @chriseclectic , it would be nice to have some additional documentation for this, e.g. a userguide or tutorial.

I think probably a userguide entry is the right format, something answering the question(s): "When and how should I use sparse representations?", and in particular how to choose the simulation frame when simulating a sparse model.

One thing to consider is whether or not we should gear this tutorial towards JAX sparse explicitly. I think yes as the scipy sparse mode still needs to be performance optimized, though the lessons/interface points hold true in this case as well.

Enact whatever actions are listed in the deprecation warnings in the 0.3.0 release notes.

Diffrax is a JAX ODE solving library providing access to many different solvers, and more control over how computations (like back propagation) are done. It also allegedly offers speed improvements over JAX's odeint function for the same method and tolerances.

What is the expected behavior?

Add access to solvers from this library through solve_ode.