Implement DynamicsBackend.solve function (#292)

Co-authored-by: DanPuzzuoli <dan.puzzuoli@gmail.com>

qiskit_dynamics/backend/dynamics_backend.py

@@ -18,13 +18,14 @@
 """
 
 import datetime
+import inspect
 import uuid
 import warnings
 
 from typing import List, Optional, Union, Dict, Tuple
 import copy
 import numpy as np
-from scipy.integrate._ivp.ivp import OdeResult  # pylint: disable=unused-import
+from scipy.integrate._ivp.ivp import OdeResult
 
 from qiskit import pulse
 from qiskit.qobj.utils import MeasLevel, MeasReturnType
@@ -43,6 +44,9 @@
 from qiskit import QiskitError, QuantumCircuit
 from qiskit import schedule as build_schedule
 from qiskit.quantum_info import Statevector, DensityMatrix
+from qiskit.quantum_info.operators.base_operator import BaseOperator
+from qiskit.quantum_info.states.quantum_state import QuantumState
+
 
 from qiskit_dynamics import RotatingFrame
 from qiskit_dynamics.array import Array
@@ -338,6 +342,56 @@ def _set_solver(self, solver):
         self._dressed_states = dressed_states
         self._dressed_states_adjoint = self._dressed_states.conj().transpose()
 
+    def solve(
+        self,
+        solve_input: List[Union[QuantumCircuit, Schedule, ScheduleBlock]],
+        t_span: Array,
+        y0: Optional[Union[Array, QuantumState, BaseOperator]] = None,
+        convert_results: Optional[bool] = True,
+        validate: Optional[bool] = True,
+    ) -> Union[OdeResult, List[OdeResult]]:
+        """Simulate a list of :class:`~qiskit.circuit.QuantumCircuit`,
+        :class:`~qiskit.pulse.Schedule`, or :class:`~qiskit.pulse.ScheduleBlock` instances and
+        return the ``OdeResult``.
+
+        This method is analogous to :meth:`.Solver.solve`, however it additionally utilizes
+        transpilation and the backend configuration to convert
+        :class:`~qiskit.circuit.QuantumCircuit` instances into pulse-level schedules for simulation.
+        The options for the solver will be drawn from ``self.options.solver_options``, and if
+        ``y0`` is not specified, it will be set from ``self.options.initial_state``.
+
+        Args:
+            t_span: Time interval to integrate over.
+            y0: Initial state.
+            solve_input: Time evolution of the system in terms of quantum circuits or qiskit
+                         pulse schedules.
+            convert_results: If ``True``, convert returned solver state results to the same class
+                             as y0. If ``False``, states will be returned in the native array type
+                             used by the specified solver method.
+            validate: Whether or not to run validation checks on the input.
+
+        Returns:
+            OdeResult: object with formatted output types.
+        """
+        if validate:
+            _validate_run_input(solve_input)
+        schedules, _ = _to_schedule_list(solve_input, backend=self)
+
+        # use default y0 if not given as parameter
+        if y0 is None:
+            y0 = self.options.initial_state
+        if y0 == "ground_state":
+            y0 = Statevector(self._dressed_states[:, 0])
+
+        solver_results = self.options.solver.solve(
+            t_span=t_span,
+            y0=y0,
+            signals=schedules,
+            convert_results=convert_results,
+            **self.options.solver_options,
+        )
+        return solver_results
+
     # pylint: disable=arguments-differ
     def run(
         self,
@@ -871,16 +925,20 @@ def default_experiment_result_function(
         raise QiskitError(f"meas_level=={backend.options.meas_level} not implemented.")
 
 
-def _validate_run_input(run_input, accept_list=True):
+def _validate_run_input(run_input, accept_list=True, caller_func_name: str = None):
     """Raise errors if the run_input is not one of QuantumCircuit, Schedule, ScheduleBlock, or
     a list of these.
     """
+    if caller_func_name is None:
+        caller_func_name = inspect.stack()[1].function
     if isinstance(run_input, list) and accept_list:
         # if list apply recursively, but no longer accept lists
         for x in run_input:
-            _validate_run_input(x, accept_list=False)
+            _validate_run_input(x, accept_list=False, caller_func_name=caller_func_name)
     elif not isinstance(run_input, (QuantumCircuit, Schedule, ScheduleBlock)):
-        raise QiskitError(f"Input type {type(run_input)} not supported by DynamicsBackend.run.")
+        raise QiskitError(
+            f"Input type {type(run_input)} not supported by DynamicsBackend.{caller_func_name}."
+        )
 
 
 def _get_acquire_instruction_timings(

releasenotes/notes/dynamics-backend-solve-fd42dac884ff80cc.yaml

@@ -0,0 +1,5 @@
+---
+features:
+  - |
+    Adds the :meth:`.DynamicsBackend.solve` method for running simulations of circuits and schedules
+    for arbitrary input types, and returning the ODE simulation results.

test/dynamics/backend/test_dynamics_backend.py

@@ -16,15 +16,17 @@
 """
 
 from types import SimpleNamespace
+from itertools import product
 
 import numpy as np
 from scipy.integrate._ivp.ivp import OdeResult
 from scipy.sparse import csr_matrix
+from scipy.linalg import expm
 
 from qiskit import QiskitError, pulse, QuantumCircuit, QuantumRegister, ClassicalRegister
-from qiskit.circuit.library import XGate, Measure
+from qiskit.circuit.library import XGate, UnitaryGate, Measure
 from qiskit.transpiler import Target, InstructionProperties
-from qiskit.quantum_info import Statevector, DensityMatrix
+from qiskit.quantum_info import Statevector, DensityMatrix, Operator, SuperOp
 from qiskit.result.models import ExperimentResult, ExperimentResultData
 from qiskit.providers.models.backendconfiguration import UchannelLO
 from qiskit.providers.backend import QubitProperties
@@ -301,6 +303,61 @@ def test_pi_pulse(self):
         counts = self.iq_to_counts(result.get_memory())
         self.assertDictEqual(counts, {"1": 1024})
 
+    def test_solve(self):
+        """Test the ODE simulation with different input and y0 types using a X pulse."""
+
+        # build solver to use in the test
+        static_ham = 2 * np.pi * 5 * np.array([[-1.0, 0.0], [0.0, 1.0]]) / 2
+        drive_op = 2 * np.pi * 0.1 * np.array([[0.0, 1.0], [1.0, 0.0]]) / 2
+
+        solver = Solver(
+            static_hamiltonian=static_ham,
+            hamiltonian_operators=[drive_op],
+            hamiltonian_channels=["d0"],
+            channel_carrier_freqs={"d0": 5.0},
+            dt=0.1,
+            rotating_frame=static_ham,
+            rwa_cutoff_freq=5.0,
+        )
+
+        backend = DynamicsBackend(solver=solver, solver_options={"atol": 1e-10, "rtol": 1e-10})
+
+        # create the circuit, pulse schedule and calibrate the gate
+        x_circ0 = QuantumCircuit(1)
+        x_circ0.x(0)
+        n_samples = 5
+        with pulse.build() as x_sched0:
+            pulse.play(pulse.Waveform([1.0] * n_samples), pulse.DriveChannel(0))
+        x_circ0.add_calibration("x", [0], x_sched0)
+
+        # create the initial states and expected simulation results
+        generator = np.array([[0, 1], [1, 0]], dtype=np.complex128)
+        rotation_strength = n_samples / 100
+        expected_unitary = expm(-1.0j * 0.5 * np.pi * rotation_strength * generator)
+        y0_and_expected_results = []
+        for y0_type in [Statevector, Operator, DensityMatrix, SuperOp]:
+            y0 = y0_type(QuantumCircuit(1))
+            expected_result = y0_type(UnitaryGate(expected_unitary))
+            y0_and_expected_results.append((y0, expected_result))
+        # y0=None defaults to Statevector
+        y0_and_expected_results.append(
+            (None, Statevector(QuantumCircuit(1)).evolve(expected_unitary))
+        )
+        input_variety = [x_sched0, x_circ0]
+
+        # solve for all combinations of input types and initial states
+        for solve_input, (y0, expected_result) in product(input_variety, y0_and_expected_results):
+            solver_results = backend.solve(
+                t_span=[0, n_samples * backend.dt],
+                y0=y0,
+                solve_input=[solve_input],
+            )
+
+            # results are always a list
+            for solver_result in solver_results:
+                self.assertTrue(solver_result.success)
+                self.assertAllClose(solver_result.y[-1], expected_result, atol=1e-8, rtol=1e-8)
+
     def test_pi_pulse_initial_state(self):
         """Test simulation of a pi pulse with a different initial state."""