Merge 2c600fe into 85303bb

doc/pulse-paper/qft.py

@@ -11,7 +11,7 @@
 
 
 import numpy as np
-from qutip import basis, fidelity
+from qutip import basis, fidelity, Options
 from qutip_qip.device import LinearSpinChain
 from qutip_qip.algorithms import qft_gate_sequence
 
@@ -23,7 +23,11 @@
 # Pulse-level simulation
 processor = LinearSpinChain(num_qubits)
 processor.load_circuit(qc)
-state2 = processor.run_state(basis([2]*num_qubits, [0]*num_qubits)).states[-1]
+options = Options(max_step=5000)
+state2 = processor.run_state(
+    basis([2]*num_qubits, [0]*num_qubits),
+    options=options
+    ).states[-1]
 
 assert(abs(1 - fidelity(state1, state2)) < 1.e-4)
 

src/qutip_qip/device/processor.py

@@ -1188,14 +1188,33 @@ def run_state(
         else:
             kwargs["c_ops"] = sys_c_ops
 
-        # choose solver:
+        # set tlist
         if "tlist" in kwargs:
             tlist = kwargs["tlist"]
             del kwargs["tlist"]
         else:
             # TODO, this can be simplified further, tlist in the solver only
             # determines the time step for intermediate result.
             tlist = self.get_full_tlist()
+        # Set the max step size as 1/10 of the total circuit time.
+        # A better solution is to use the gate, which
+        # is however, much harder to implement at this stage, see also
+        # https://github.com/qutip/qutip-qip/issues/184.
+        full_tlist = self.get_full_tlist()
+        if full_tlist is not None:
+            total_circuit_time = (full_tlist)[-1]
+        else:
+            total_circuit_time = 0.
+        if is_qutip5:
+            options = kwargs.get("options", qutip.SolverOptions())
+            if options["max_step"] == 0.0:
+                options["max_step"] = total_circuit_time / 10
+        else:
+            options = kwargs.get("options", qutip.Options())
+            if options.max_step == 0.0:
+                options.max_step = total_circuit_time / 10
+        kwargs["options"] = options
+        # choose solver:
         if solver == "mesolve":
             evo_result = mesolve(
                 H=noisy_qobjevo, rho0=init_state, tlist=tlist, **kwargs

tests/test_processor.py

@@ -11,7 +11,7 @@
     from qutip import Options as SolverOptions
 else:
     from qutip import SolverOptions
-from qutip_qip.device.processor import Processor
+from qutip_qip.device import Processor, LinearSpinChain
 from qutip import (
     basis, sigmaz, sigmax, sigmay, identity, destroy, tensor,
     rand_ket, rand_dm, fidelity)
@@ -20,6 +20,11 @@
     DecoherenceNoise, RandomNoise, ControlAmpNoise)
 from qutip_qip.qubits import qubit_states
 from qutip_qip.pulse import Pulse
+from qutip_qip.circuit import QubitCircuit
+if parse_version(qutip.__version__) < parse_version('5.dev'):
+    from qutip import Options as SolverOptions
+else:
+    from qutip import SolverOptions
 
 
 class TestCircuitProcessor:
@@ -383,3 +388,31 @@ def test_pulse_mode(self):
         processor.pulse_mode = "discrete"
         assert(processor.pulse_mode == "discrete")
         assert(processor.pulses[0].spline_kind == "step_func")
+
+    def test_max_step_size(self):
+        num_qubits = 2
+        init_state = tensor([basis(2, 1), basis(2, 1)])
+        qc = QubitCircuit(2)
+
+        # ISWAP acts trivially on the initial states.
+        # If no max_step are defined,
+        # the solver will choose a step size too large
+        # such that the X gate will be skipped.
+        qc.add_gate("ISWAP", targets=[0, 1])
+        qc.add_gate("ISWAP", targets=[0, 1])
+        qc.add_gate("X", targets=[0])
+        processor = LinearSpinChain(num_qubits)
+        processor.load_circuit(qc)
+
+        # No max_step
+        final_state = processor.run_state(
+            init_state,
+            options=SolverOptions(max_step=10000)  # too large max_step
+        ).states[-1]
+        expected_state = tensor([basis(2, 0), basis(2, 1)])
+        assert pytest.approx(fidelity(final_state, expected_state), 0.001) == 0
+
+        # With default max_step
+        final_state = processor.run_state(init_state).states[-1]
+        expected_state = tensor([basis(2, 0), basis(2, 1)])
+        assert pytest.approx(fidelity(final_state, expected_state), 0.001) == 1