Closing in

Author: pmenczel

qutip/solver/integrator/scipy_integrator.py

@@ -298,15 +298,10 @@ def set_state(self, t, state0):
         self._size = state0.shape[0]
         if self._mat_state:
             state0 = _data.column_stack(state0)
-        print(' --- HERE --- ')
-        print(state0)
-        print(state0.to_array().ravel().view(np.float64))
         self._ode_solver.set_initial_value(
             state0.to_array().ravel().view(np.float64),
             t
         )
-        print(self._ode_solver._y)
-        print(' --- --- ')
 
     def _check_failed_integration(self):
         if self._ode_solver.successful():
@@ -413,6 +408,7 @@ def mcstep(self, t, copy=True):
         # a range in which the state at any time can be computed with
         # minimal work. We keep track of the range with _back[0] and _front.
         self._check_handle()
+        print("mcstep: t=", t)
         if self._ode_solver.t == t:
             # Exact same `t` as the last call, nothing to do.
             pass
@@ -444,20 +440,29 @@ def _one_step(self, t):
         # `rhs` can cause exceptions to this, but _backstep catch those cases.
         safe_delta = self._ode_solver._integrator.rwork[11]/100 + 1e-15
         t_ode = self._ode_solver.t
+        if t_ode > 0.16:
+            print('one_step: t=', t ,', t_ode=', t_ode)
 
         if t > self._front and t_ode >= self._front:
             # The state is at self._front, do a step
             self._back = self.get_state()
-            print(' --- THERE --- ', t_ode, ' ,  ', t)
-            print(self._ode_solver._y)
+            if t_ode > 0.16:
+                print(' - case 1: ', np.isnan(self._ode_solver._y).any())
             self._ode_solver.integrate(min(self._front + safe_delta, t))
-            print(' --- --- ')
+            if t_ode > 0.16:
+                print(' - after integrate (', self._ode_solver.get_return_code(), '): ', np.isnan(self._ode_solver._y).any())
             self._front = self._ode_solver._integrator.rwork[12]
             # We asked for a fraction of a step, now complete it.
             self._ode_solver.integrate(min(self._front, t))
+            if t_ode > 0.16:
+                print(' - after integrate 2 (', self._ode_solver.get_return_code(), '): ', np.isnan(self._ode_solver._y).any())
         elif t > self._front:
             # The state is at a time before t_front, advance to t_front
+            if t_ode > 0.16:
+                print(' - case 2: ', np.isnan(self._ode_solver._y).any())
             self._ode_solver.integrate(self._front)
+            if t_ode > 0.16:
+                print(' - after integrate (', self._ode_solver.get_return_code(), '): ', np.isnan(self._ode_solver._y).any())
 
     def _backstep(self, t):
         """

qutip/solver/multitraj.py

@@ -138,9 +138,7 @@ def step(
         """
         if not self._integrator._is_set:
             raise RuntimeError("The `start` method must called first.")
-        print('before args change: ', self._integrator._integrator._ode_solver._y)
         self._argument(args)
-        print('after args change: ', self._integrator._integrator._ode_solver._y)
         _, state = self._integrator.integrate(t, copy=False)
         return self._restore_state(state, copy=copy)
 

qutip/tests/solver/test_mcsolve.py

@@ -516,14 +516,11 @@ def test_MCSolver_stepping():
     H = qutip.num(size)
     solver = MCSolver(H, a)
     solver.start(qutip.basis(size, size-1), 0, seed=0)
-    print('a', solver._integrator._integrator._ode_solver._y)
     solver.options = {'method': 'lsoda'}
-    print('b', solver._integrator._integrator._ode_solver._y)
     state = solver.step(1)
     assert qutip.expect(qutip.qeye(size), state) == pytest.approx(1)
     assert qutip.expect(qutip.num(size), state) == pytest.approx(size - 1)
     assert state.isket
-    print('c', solver._integrator._integrator._ode_solver._y)
     state = solver.step(5, args={'coupling': 5})
     assert qutip.expect(qutip.qeye(size), state) == pytest.approx(1)
     assert qutip.expect(qutip.num(size), state) <= size - 1