Merge branch 'master' into feature/clean-up-bloch-sphere-rendering

.github/workflows/tests.yml

@@ -27,6 +27,7 @@ jobs:
         # matrix size; make sure to test all supported versions in some form.
         python-version: [3.9]
         case-name: [defaults]
+        numpy-requirement: [">=1.20,<1.21"]
         # Extra special cases.  In these, the new variable defined should always
         # be a truth-y value (hence 'nomkl: 1' rather than 'mkl: 0'), because
         # the lack of a variable is _always_ false-y, and the defaults lack all
@@ -38,6 +39,8 @@ jobs:
             os: ubuntu-latest
             python-version: 3.6
             scipy-requirement: ">=1.4,<1.5"
+            # let the old SciPy version select an appropriate numpy version:
+            numpy-requirement: ""
 
           # No MKL runs.  MKL is now the default for conda installations, but
           # not necessarily for pip.
@@ -51,6 +54,7 @@ jobs:
           - case-name: OpenMP
             os: ubuntu-latest
             python-version: 3.9
+            numpy-requirement: ">=1.20,<1.21"
             openmp: 1
 
           # Builds without Cython at runtime.  This is a core feature;
@@ -79,13 +83,13 @@ jobs:
           fi
           export CI_QUTIP_WITH_OPENMP=${{ matrix.openmp }}
           if [[ -z "${{ matrix.nomkl }}" ]]; then
-            conda install blas=*=mkl numpy "scipy${{ matrix.scipy-requirement }}"
+            conda install blas=*=mkl "numpy${{ matrix.numpy-requirement }}" "scipy${{ matrix.scipy-requirement }}"
           elif [[ "${{ matrix.os }}" =~ ^windows.*$ ]]; then
             # Conda doesn't supply forced nomkl builds on Windows, so we rely on
             # pip not automatically linking to MKL.
-            pip install numpy "scipy${{ matrix.scipy-requirement }}"
+            pip install "numpy${{ matrix.numpy-requirement }}" "scipy${{ matrix.scipy-requirement }}"
           else
-            conda install nomkl numpy "scipy${{ matrix.scipy-requirement }}"
+            conda install nomkl "numpy${{ matrix.numpy-requirement }}" "scipy${{ matrix.scipy-requirement }}"
           fi
           python -m pip install -e .[$QUTIP_TARGET]
           python -m pip install pytest-cov coveralls
@@ -95,6 +99,15 @@ jobs:
           conda list
           python -c "import qutip; qutip.about()"
 
+      - name: Environment information
+        run: |
+          uname -a
+          if [[ "ubuntu-latest" == "${{ matrix.os }}" ]]; then
+              hostnamectl
+              lscpu
+              free -h
+          fi
+
       - name: Run tests
         # If our tests are running for longer than an hour, _something_ is wrong
         # somewhere.  The GitHub default is 6 hours, which is a bit long to wait

doc/apidoc/functions.rst

@@ -334,7 +334,7 @@ Utility Functions
 -----------------
 
 .. automodule:: qutip.utilities
-    :members: n_thermal, linspace_with, clebsch, convert_unit
+    :members: n_thermal, clebsch, convert_unit
 
 
 .. _functions-fileio:

doc/conf.py

@@ -227,6 +227,10 @@ def qutip_version():
 # This is the file name suffix for HTML files (e.g. ".xhtml").
 #html_file_suffix = None
 
+html_css_files = [
+    'site.css',
+]
+
 # -- Options for HTMLHelp output ------------------------------------------
 
 # Output file base name for HTML help builder.

doc/contributors.rst

@@ -79,7 +79,7 @@ Contributors
 - Louis Tessler
 - Lucas Verney
 - Marco David
-- Marek
+- Marek Narozniak
 - Markus Baden
 - Martín Sande
 - Mateo Laguna

doc/guide/dynamics/dynamics-floquet.rst

@@ -206,7 +206,7 @@ For convenience, all the steps described above for calculating the evolution of
 
 .. code-block:: python
 
-    output = fsesolve(H, psi0, times, [num(2)], args)
+    output = fsesolve(H, psi0=psi0, tlist=tlist, e_ops=[qutip.num(2)], args=args)
     p_ex = output.expect[0]
 
 .. _floquet-dissipative:

doc/guide/dynamics/dynamics-monte.rst

@@ -39,16 +39,18 @@ If more than a single collapse operator is present in Eq. :eq:`heff`, the probab
 
 Evaluating the MC evolution to first-order in time is quite tedious.  Instead, QuTiP uses the following algorithm to simulate a single realization of a quantum system.  Starting from a pure state :math:`\left|\psi(0)\right>`:
 
-- **I:** Choose a random number :math:`r` between zero and one, representing the probability that a quantum jump occurs.
+- **Ia:** Choose a random number :math:`r_1` between zero and one, representing the probability that a quantum jump occurs.
 
-- **II:** Integrate the Schrödinger equation, using the effective Hamiltonian :eq:`heff` until a time :math:`\tau` such that the norm of the wave function satisfies :math:`\left<\psi(\tau)\right.\left|\psi(\tau)\right>=r`, at which point a jump occurs.
+- **Ib:** Choose a random number :math:`r_2` between zero and one, used to select which collapse operator was responsible for the jump.
+
+- **II:** Integrate the Schrödinger equation, using the effective Hamiltonian :eq:`heff` until a time :math:`\tau` such that the norm of the wave function satisfies :math:`\left<\psi(\tau)\right.\left|\psi(\tau)\right> = r_1`, at which point a jump occurs.
 
 - **III:** The resultant jump projects the system at time :math:`\tau` into one of the renormalized states given by Eq. :eq:`project`.  The corresponding collapse operator :math:`C_{n}` is chosen such that :math:`n` is the smallest integer satisfying:
 
 .. math::
     :label: mc3
 
-    \sum_{i=1}^{n} P_{n}(\tau) \ge r
+    \sum_{i=1}^{n} P_{n}(\tau) \ge r_2
 
 where the individual :math:`P_{n}` are given by Eq. :eq:`pcn`.  Note that the left hand side of Eq. :eq:`mc3` is, by definition, normalized to unity.
 

doc/requirements.txt

@@ -1,6 +1,6 @@
 alabaster==0.7.12
 appnope==0.1.2
-Babel==2.9.0
+Babel==2.9.1
 backcall==0.2.0
 certifi==2020.12.5
 chardet==4.0.0
@@ -23,7 +23,7 @@ packaging==20.9
 parso==0.8.2
 pexpect==4.8.0
 pickleshare==0.7.5
-Pillow==8.2.0
+Pillow==8.3.2
 prompt-toolkit==3.0.18
 ptyprocess==0.7.0
 Pygments==2.8.1

doc/static/site.css

@@ -1,64 +1,19 @@
-@import url('https://fonts.googleapis.com/css?family=Source+Code+Pro');
-
-
-.navbar-text {
-    color: #e8e8e8 !important;
-}
-
-a {
-  color: #599AD3;
-  text-decoration: none;
-}
-
-a:hover,
-a:focus {
-  color: #8C0028;
-  text-decoration: underline;
-}
-
-a:focus {
-  outline: thin dotted #333;
-  outline: 5px auto -webkit-focus-ring-color;
-  outline-offset: -2px;
-}
-
-code,
-pre {
-  padding: 0 3px 2px;
-  font-family: "Source Code Pro", Monaco, Menlo, Consolas, "Courier New", monospace;
-  font-size: 12px;
-  color: #333333;
-  -webkit-border-radius: 3px;
-     -moz-border-radius: 3px;
-          border-radius: 3px;
-}
-
-.alert {
-  border-width: 2px;
-  color: #09224F;
-  font-weight: bold;
-  -webkit-border-radius: 3px;
-  -moz-border-radius: 3px;
-       border-radius: 3px;
-  }
-
-.alert-success {
-  background-color: #B9E0B0;
-  border-color: #79C36A;
-}
-
-
-.alert-info {
-  background-color: #A6CBE9;
-  border-color: #599AD3;
+/* Fix for: https://github.com/readthedocs/sphinx_rtd_theme/issues/301 */
+/* Fix taken from: https://github.com/readthedocs/sphinx_rtd_theme/pull/383/ */
+span.eqno {
+  margin-left: 5px;
+  float: right;
+  /* position the number above the equation so that :hover is activated */
+  z-index: 1;
+  position: relative;
 }
 
-.alert-warning {
- background-color: #FBD1A7;
-   border-color: #F9A65A;
+span.eqno .headerlink {
+  display: none;
+  visibility: hidden;
 }
 
-.alert-danger {
-background-color: #F7A7AA;
-  border-color: #F1595F;
+span.eqno:hover .headerlink {
+  display: inline-block;
+  visibility: visible;
 }

qutip/control/dynamics.py

@@ -185,9 +185,7 @@ class Dynamics(object):
 
     oper_dtype : type
         Data type for internal dynamics generators, propagators and time
-        evolution operators. This can be ndarray or Qobj, or (in theory) any
-        other representaion that supports typical matrix methods (e.g. dot)
-        ndarray performs best for smaller quantum systems.
+        evolution operators. This can be ndarray or Qobj.
         Qobj may perform better for larger systems, and will also
         perform better when (custom) fidelity measures use Qobj methods
         such as partial trace.
@@ -386,7 +384,7 @@ def reset(self):
         self.time_depend_ctrl_dyn_gen = False
         # These internal attributes will be of the internal operator data type
         # used to compute the evolution
-        # Note this maybe ndarray, Qobj or some other depending on oper_dtype
+        # This will be either ndarray or Qobj
         self._drift_dyn_gen = None
         self._ctrl_dyn_gen = None
         self._phased_ctrl_dyn_gen = None
@@ -684,7 +682,7 @@ def _choose_oper_dtype(self):
             else:
                 ctrls = self.ctrl_dyn_gen
             for c in ctrls:
-               dg = dg + c
+                dg = dg + c
 
             N = dg.data.shape[0]
             n = dg.data.nnz
@@ -724,7 +722,15 @@ def _init_evo(self):
         self.sys_shape = self.initial.shape
         # Set the phase application method
         self._init_phase()
+
         self._set_memory_optimizations()
+        if self.sparse_eigen_decomp and self.sys_shape[0] <= 2:
+            raise ValueError(
+                "Single qubit pulse optimization dynamics cannot use sparse"
+                " eigenvector decomposition because of limitations in"
+                " scipy.linalg.eigsh. Pleae set sparse_eigen_decomp to False"
+                " or increase the size of the system.")
+
         n_ts = self.num_tslots
         n_ctrls = self.num_ctrls
         if self.oper_dtype == Qobj:
@@ -748,26 +754,10 @@ def _init_evo(self):
             else:
                 self._ctrl_dyn_gen = [ctrl.full()
                                         for ctrl in self.ctrl_dyn_gen]
-        elif self.oper_dtype == sp.csr_matrix:
-            self._initial = self.initial.data
-            self._target = self.target.data
-            if self.time_depend_drift:
-                self._drift_dyn_gen = [d.data for d in self.drift_dyn_gen]
-            else:
-                self._drift_dyn_gen = self.drift_dyn_gen.data
-
-            if self.time_depend_ctrl_dyn_gen:
-                self._ctrl_dyn_gen = np.empty([n_ts, n_ctrls], dtype=object)
-                for k in range(n_ts):
-                    for j in range(n_ctrls):
-                        self._ctrl_dyn_gen[k, j] = \
-                                    self.ctrl_dyn_gen[k, j].data
-            else:
-                self._ctrl_dyn_gen = [ctrl.data for ctrl in self.ctrl_dyn_gen]
         else:
-            logger.warn("Unknown option '{}' for oper_dtype. "
-                "Assuming that internal drift, ctrls, initial and target "
-                "have been set correctly".format(self.oper_dtype))
+            raise ValueError(
+                "Unknown oper_dtype {!r}. The oper_dtype may be qutip.Qobj or"
+                " numpy.ndarray.".format(self.oper_dtype))
 
         if self.cache_phased_dyn_gen:
             if self.time_depend_ctrl_dyn_gen:
@@ -1149,24 +1139,19 @@ def _get_onto_evo_target(self):
             #Target is operator
             targ = la.inv(self.target.full())
             if self.oper_dtype == Qobj:
-                self._onto_evo_target = Qobj(targ)
+                rev_dims = [self.target.dims[1], self.target.dims[0]]
+                self._onto_evo_target = Qobj(targ, dims=rev_dims)
             elif self.oper_dtype == np.ndarray:
                 self._onto_evo_target = targ
-            elif self.oper_dtype == sp.csr_matrix:
-                self._onto_evo_target = sp.csr_matrix(targ)
             else:
-                targ_cls = self._target.__class__
-                self._onto_evo_target = targ_cls(targ)
+                assert False, f"Unknown oper_dtype {self.oper_dtype!r}"
         else:
             if self.oper_dtype == Qobj:
                 self._onto_evo_target = self.target.dag()
             elif self.oper_dtype == np.ndarray:
                 self._onto_evo_target = self.target.dag().full()
-            elif self.oper_dtype == sp.csr_matrix:
-                self._onto_evo_target = self.target.dag().data
             else:
-                targ_cls = self._target.__class__
-                self._onto_evo_target = targ_cls(self.target.dag().full())
+                assert False, f"Unknown oper_dtype {self.oper_dtype!r}"
 
         return self._onto_evo_target
 
@@ -1570,18 +1555,19 @@ def _spectral_decomp(self, k):
             eig_val, eig_vec = sp_eigs(H.data, H.isherm,
                                        sparse=self.sparse_eigen_decomp)
             eig_vec = eig_vec.T
+            if self.sparse_eigen_decomp:
+                # when sparse=True, sp_eigs returns an ndarray where each
+                # element is a sparse matrix so we convert it into a sparse
+                # matrix we can later pass to Qobj(...)
+                eig_vec = sp.hstack(eig_vec)
 
         elif self.oper_dtype == np.ndarray:
             H = self._dyn_gen[k]
             # returns row vector of eigenvals, columns with the eigenvecs
             eig_val, eig_vec = eigh(H)
+
         else:
-            if sparse:
-                H = self._dyn_gen[k].toarray()
-            else:
-                H = self._dyn_gen[k]
-            # returns row vector of eigenvals, columns with the eigenvecs
-            eig_val, eig_vec = eigh(H)
+            assert False, f"Unknown oper_dtype {self.oper_dtype!r}"
 
         # assuming H is an nxn matrix, find n
         n = self.get_drift_dim()
@@ -1630,14 +1616,16 @@ def _spectral_decomp(self, k):
             if self._dyn_gen_eigenvectors_adj is not None:
                 self._dyn_gen_eigenvectors_adj[k] = \
                             self._dyn_gen_eigenvectors[k].dag()
-        else:
+        elif self.oper_dtype == np.ndarray:
             self._prop_eigen[k] = np.diagflat(prop_eig)
             self._dyn_gen_eigenvectors[k] = eig_vec
             # The _dyn_gen_eigenvectors_adj list is not used in
             # memory optimised modes
             if self._dyn_gen_eigenvectors_adj is not None:
                 self._dyn_gen_eigenvectors_adj[k] = \
                             self._dyn_gen_eigenvectors[k].conj().T
+        else:
+            assert False, f"Unknown oper_dtype {self.oper_dtype!r}"
 
     def _get_dyn_gen_eigenvectors_adj(self, k):
         # The _dyn_gen_eigenvectors_adj list is not used in
@@ -1736,10 +1724,8 @@ def _get_omega(self):
             if self.oper_dtype == Qobj:
                 self._omega = Qobj(omg, dims=self.dyn_dims)
                 self._omega_qobj = self._omega
-            elif self.oper_dtype == sp.csr_matrix:
-                self._omega = sp.csr_matrix(omg)
             else:
-                 self._omega = omg
+                self._omega = omg
         return self._omega
 
     def _set_phase_application(self, value):