Merge pull request #183 from BoxiLi/max_step

Add a default max_step based on the total circuit time.

This does not solve the problem completely. But it will avoid those mistakes that keep popping up when people try and compile one gate and send it to the solver. When there are more gates and a non-trivial circuit, the probability of this happening will be very low.

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.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