Merge pull request #1565 from quantshah/prepare-4.6.2

Prepare 4.6.2 release

README.md

@@ -76,3 +76,9 @@ You can also help out with users' questions, or discuss proposed changes in the
 All code contributions are acknowledged in the [contributors](http://qutip.org/docs/latest/contributors.html) section in the documentation.
 
 For more information, including technical advice, please see the ["contributing to QuTiP development" section of the documentation](http://qutip.org/docs/latest/development/contributing.html).
+
+
+Citing QuTiP
+------------
+
+If you use QuTiP in your research, please cite the original QuTiP papers that are available [here](https://dml.riken.jp/?s=QuTiP).

VERSION

@@ -1 +1 @@
-4.6.1
+4.6.2

doc/apidoc/functions.rst

@@ -99,7 +99,7 @@ Entropy Functions
 -----------------
 
 .. automodule:: qutip.entropy
-    :members: concurrence, entropy_conditional, entropy_linear, entropy_mutual, entropy_vn
+    :members: concurrence, entropy_conditional, entropy_linear, entropy_mutual, entropy_relative, entropy_vn
 
 
 Density Matrix Metrics

doc/changelog.rst

@@ -7,6 +7,33 @@
 Change Log
 **********
 
+Version 4.6.2 (June 2, 2021)
+++++++++++++++++++++++++++++
+
+This minor release adds a function to calculate the quantum relative entropy, fixes a corner case in handling time-dependent Hamiltonians in ``mesolve`` and adds back support for a wider range of matplotlib versions when plotting or animating Bloch spheres.
+
+It also adds a section in the README listing the papers which should be referenced while citing QuTiP. 
+
+
+Improvements
+------------
+- Added a "Citing QuTiP" section to the README, containing a link to the QuTiP papers. (`#1554 <https://github.com/qutip/qutip/pull/1554>`_)
+- Added ``entropy_relative`` which returns the quantum relative entropy between two density matrices. (`#1553 <https://github.com/qutip/qutip/pull/1553>`_)
+
+Bug Fixes
+---------
+- Fixed Bloch sphere distortion when using Matplotlib >= 3.3.0. (`#1496  <https://github.com/qutip/qutip/pull/1496>`_)
+- Removed use of integer-like floats in math.factorial since it is deprecated as of Python 3.9. (`#1550 <https://github.com/qutip/qutip/pull/1550>`_)
+- Simplified call to ffmpeg used in the the Bloch sphere animation tutorial to work with recent versions of ffmpeg. (`#1557 <https://github.com/qutip/qutip/pull/1557>`_)
+- Removed blitting in Bloch sphere FuncAnimation example. (`#1558 <https://github.com/qutip/qutip/pull/1558>`_)
+- Added a version checking condition to handle specific functionalities depending on the matplotlib version. (`#1556 <https://github.com/qutip/qutip/pull/1556>`_)
+- Fixed ``mesolve`` handling of time-dependent Hamiltonian with a custom tlist and ``c_ops``. (`#1561 <https://github.com/qutip/qutip/pull/1561>`_)
+
+Developer Changes
+-----------------
+- Read documentation version and release from the VERSION file.
+
+
 Version 4.6.1 (May 4, 2021)
 +++++++++++++++++++++++++++
 

doc/conf.py

@@ -5,8 +5,9 @@
 # add these directories to sys.path here. If the directory is relative to the
 # documentation root, use os.path.abspath to make it absolute, like shown here.
 #
-import sys, os
-sys.path.append(os.path.abspath('sphinxext'))
+import os
+import pathlib
+import warnings
 
 # -- General configuration ------------------------------------------------
 
@@ -65,14 +66,51 @@
 
 copyright = '2011 and later, {}'.format(author)
 
+
+def _check_source_folder_and_imported_qutip_match():
+    """ Warn if the imported qutip and the source folder the documentation
+        is being built from don't match.
+
+        The generated documentation contains material from both the
+        source folder (e.g. ``.rst`` files) and from the imported qutip
+        (e.g. docstrings), so if the two don't match the generated
+        documentation will be a chimera.
+    """
+    import qutip
+    qutip_folder = pathlib.Path(qutip.__file__).absolute().parent.parent
+    source_folder = pathlib.Path(__file__).absolute().parent.parent
+    if qutip_folder != source_folder:
+        warnings.warn(
+            "The documentation source and imported qutip package are"
+            " not from the same source folder. This may result in the"
+            " documentation containing text from different sources."
+            " Documentation source: {!r}."
+            " Qutip package source: {!r}.".format(
+                str(source_folder), str(qutip_folder)
+            )
+        )
+
+
+_check_source_folder_and_imported_qutip_match()
+
+
+def qutip_version():
+    """ Retrieve the QuTiP version from ``../VERSION``.
+    """
+    src_folder_root = pathlib.Path(__file__).absolute().parent.parent
+    version = src_folder_root.joinpath(
+        "VERSION"
+    ).read_text().strip()
+    return version
+
+
 # The version info for the project you're documenting, acts as replacement for
 # |version| and |release|, also used in various other places throughout the
 # built documents.
-#
-# The short X.Y version.
-version = '4.7'
 # The full version, including alpha/beta/rc tags.
-release = '4.7.0.dev0'
+release = qutip_version()
+# The short X.Y version.
+version = ".".join(release.split(".")[:2])
 # The language for content autogenerated by Sphinx. Refer to documentation
 # for a list of supported languages.
 #
@@ -201,7 +239,7 @@
 # support for the 'physics' package.  MathJax 3 does, so once Sphinx is using
 # that (should be in Sphinx 4), you will be able to swap to using that.  In the
 # meantime, we just have to define all the functions we're going to use.
-# 
+#
 # See:
 # - http://docs.mathjax.org/en/v3.0-latest/input/tex/extensions/physics.html
 mathjax_config = {

doc/frontmatter.rst

@@ -43,7 +43,7 @@ This document contains a user guide and automatically generated API documentatio
 
 :Author: Simon Cross
 
-:version: |version|
+:release: |release|
 
 :copyright: This documentation is licensed under the Creative Commons Attribution 3.0 Unported License.
 

doc/guide/guide-bloch.rst

@@ -396,7 +396,7 @@ Generating Images for Animation
 An example of generating images for generating an animation outside of Python is given below::
 
      import numpy as np
-     b = Bloch()
+     b = qutip.Bloch()
      b.vector_color = ['r']
      b.view = [-40, 30]
      for i in range(len(sx)):
@@ -407,7 +407,7 @@ An example of generating images for generating an animation outside of Python is
 
 Generating an animation using ffmpeg (for example) is fairly simple::
 
-   ffmpeg -r 20 -b 1800 -i bloch_%01d.png bloch.mp4
+   ffmpeg -i temp/bloch_%01d.png bloch.mp4
 
 .. _bloch-animate-decay-direct:
 
@@ -425,7 +425,7 @@ The code to directly generate an mp4 movie of the Qubit decay is as follows ::
    from mpl_toolkits.mplot3d import Axes3D
 
    fig = pyplot.figure()
-   ax = Axes3D(fig,azim=-40,elev=30)
+   ax = Axes3D(fig, azim=-40, elev=30)
    sphere = qutip.Bloch(axes=ax)
 
    def animate(i):
@@ -440,7 +440,7 @@ The code to directly generate an mp4 movie of the Qubit decay is as follows ::
       return ax
 
    ani = animation.FuncAnimation(fig, animate, np.arange(len(sx)),
-                                 init_func=init, blit=True, repeat=False)
+                                 init_func=init, blit=False, repeat=False)
    ani.save('bloch_sphere.mp4', fps=20)
 
 The resulting movie may be viewed here: `bloch_decay.mp4 <https://raw.githubusercontent.com/qutip/qutip/master/doc/figures/bloch_decay.mp4>`_

qutip/bloch.py

@@ -38,16 +38,29 @@
 from numpy import (ndarray, array, linspace, pi, outer, cos, sin, ones, size,
                    sqrt, real, mod, append, ceil, arange)
 
+from packaging.version import parse as parse_version
+
 from qutip.qobj import Qobj
 from qutip.expect import expect
 from qutip.operators import sigmax, sigmay, sigmaz
 
 try:
+    import matplotlib
     import matplotlib.pyplot as plt
     from mpl_toolkits.mplot3d import Axes3D
     from matplotlib.patches import FancyArrowPatch
     from mpl_toolkits.mplot3d import proj3d
 
+    # Define a custom _axes3D function based on the matplotlib version.
+    # The auto_add_to_figure keyword is new for matplotlib>=3.4.
+    if parse_version(matplotlib.__version__) >= parse_version('3.4'):
+        def _axes3D(fig, *args, **kwargs):
+            ax = Axes3D(fig, *args, auto_add_to_figure=False, **kwargs)
+            return fig.add_axes(ax)
+    else:
+        def _axes3D(*args, **kwargs):
+            return Axes3D(*args, **kwargs)
+
     class Arrow3D(FancyArrowPatch):
         def __init__(self, xs, ys, zs, *args, **kwargs):
             FancyArrowPatch.__init__(self, (0, 0), (0, 0), *args, **kwargs)
@@ -56,7 +69,7 @@ def __init__(self, xs, ys, zs, *args, **kwargs):
 
         def draw(self, renderer):
             xs3d, ys3d, zs3d = self._verts3d
-            xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, renderer.M)
+            xs, ys, zs = proj3d.proj_transform(xs3d, ys3d, zs3d, self.axes.M)
 
             self.set_positions((xs[0], ys[0]), (xs[1], ys[1]))
             FancyArrowPatch.draw(self, renderer)
@@ -448,7 +461,8 @@ def render(self, fig=None, axes=None):
             self.fig = plt.figure(figsize=self.figsize)
 
         if not axes:
-            self.axes = Axes3D(self.fig, azim=self.view[0], elev=self.view[1])
+            self.axes = _axes3D(self.fig, azim=self.view[0],
+                                elev=self.view[1])
 
         if self.background:
             self.axes.clear()
@@ -461,6 +475,10 @@ def render(self, fig=None, axes=None):
             self.axes.set_xlim3d(-0.7, 0.7)
             self.axes.set_ylim3d(-0.7, 0.7)
             self.axes.set_zlim3d(-0.7, 0.7)
+        # Manually set aspect ratio to fit a square bounding box.
+        # Matplotlib did this stretching for < 3.3.0, but not above.
+        if parse_version(matplotlib.__version__) >= parse_version('3.3'):
+            self.axes.set_box_aspect((1, 1, 1))
 
         self.axes.grid(False)
         self.plot_back()

qutip/entropy.py

@@ -32,9 +32,10 @@
 ###############################################################################
 
 __all__ = ['entropy_vn', 'entropy_linear', 'entropy_mutual', 'negativity',
-           'concurrence', 'entropy_conditional', 'entangling_power']
+           'concurrence', 'entropy_conditional', 'entangling_power',
+           'entropy_relative']
 
-from numpy import e, real, sort, sqrt
+from numpy import conj, e, inf, imag, inner, real, sort, sqrt
 from numpy.lib.scimath import log, log2
 from qutip.qobj import ptrace
 from qutip.states import ket2dm
@@ -223,45 +224,98 @@ def entropy_mutual(rho, selA, selB, base=e, sparse=False):
     return out
 
 
-def _entropy_relative(rho, sigma, base=e, sparse=False):
+def entropy_relative(rho, sigma, base=e, sparse=False, tol=1e-12):
     """
-    ****NEEDS TO BE WORKED ON****
-
     Calculates the relative entropy S(rho||sigma) between two density
     matrices.
 
     Parameters
     ----------
-    rho : qobj
-        First density matrix.
-    sigma : qobj
-        Second density matrix.
+    rho : :class:`qutip.Qobj`
+        First density matrix (or ket which will be converted to a density
+        matrix).
+    sigma : :class:`qutip.Qobj`
+        Second density matrix (or ket which will be converted to a density
+        matrix).
     base : {e,2}
-        Base of logarithm.
+        Base of logarithm. Defaults to e.
+    sparse : bool
+        Flag to use sparse solver when determining the eigenvectors
+        of the density matrices. Defaults to False.
+    tol : float
+        Tolerance to use to detect 0 eigenvalues or dot producted between
+        eigenvectors. Defaults to 1e-12.
 
     Returns
     -------
     rel_ent : float
-        Value of relative entropy.
+        Value of relative entropy. Guaranteed to be greater than zero
+        and should equal zero only when rho and sigma are identical.
+
+    Examples
+    --------
+
+    First we define two density matrices:
+
+    >>> rho = qutip.ket2dm(qutip.ket("00"))
+    >>> sigma = rho + qutip.ket2dm(qutip.ket("01"))
+    >>> sigma = sigma.unit()
 
+    Then we calculate their relative entropy using base 2 (i.e. ``log2``)
+    and base e (i.e. ``log``).
+
+    >>> qutip.entropy_relative(rho, sigma, base=2)
+    1.0
+    >>> qutip.entropy_relative(rho, sigma)
+    0.6931471805599453
+
+    References
+    ----------
+
+    See Nielsen & Chuang, "Quantum Computation and Quantum Information",
+    Section 11.3.1, pg. 511 for a detailed explanation of quantum relative
+    entropy.
     """
-    if rho.type != 'oper' or sigma.type != 'oper':
-        raise TypeError("Inputs must be density matrices..")
-    # sigma terms
-    svals = sp_eigs(sigma.data, sigma.isherm, vecs=False, sparse=sparse)
-    snzvals = svals[svals != 0]
+    if rho.isket:
+        rho = ket2dm(rho)
+    if sigma.isket:
+        sigma = ket2dm(sigma)
+    if not rho.isoper or not sigma.isoper:
+        raise TypeError("Inputs must be density matrices.")
+    if rho.dims != sigma.dims:
+        raise ValueError("Inputs must have the same shape and dims.")
     if base == 2:
-        slogvals = log2(snzvals)
+        log_base = log2
     elif base == e:
-        slogvals = log(snzvals)
+        log_base = log
     else:
         raise ValueError("Base must be 2 or e.")
-    # rho terms
-    rvals = sp_eigs(rho.data, rho.isherm, vecs=False, sparse=sparse)
-    rnzvals = rvals[rvals != 0]
-    # calculate tr(rho*log sigma)
-    rel_trace = float(real(sum(rnzvals * slogvals)))
-    return -entropy_vn(rho, base, sparse) - rel_trace
+    # S(rho || sigma) = sum_i(p_i log p_i) - sum_ij(p_i P_ij log q_i)
+    #
+    # S is +inf if the kernel of sigma (i.e. svecs[svals == 0]) has non-trivial
+    # intersection with the support of rho (i.e. rvecs[rvals != 0]).
+    rvals, rvecs = sp_eigs(rho.data, rho.isherm, vecs=True, sparse=sparse)
+    if any(abs(imag(rvals)) >= tol):
+        raise ValueError("Input rho has non-real eigenvalues.")
+    rvals = real(rvals)
+    svals, svecs = sp_eigs(sigma.data, sigma.isherm, vecs=True, sparse=sparse)
+    if any(abs(imag(svals)) >= tol):
+        raise ValueError("Input sigma has non-real eigenvalues.")
+    svals = real(svals)
+    # Calculate inner products of eigenvectors and return +inf if kernel
+    # of sigma overlaps with support of rho.
+    P = abs(inner(rvecs, conj(svecs))) ** 2
+    if (rvals >= tol) @ (P >= tol) @ (svals < tol):
+        return inf
+    # Avoid -inf from log(0) -- these terms will be multiplied by zero later
+    # anyway
+    svals[abs(svals) < tol] = 1
+    nzrvals = rvals[abs(rvals) >= tol]
+    # Calculate S
+    S = nzrvals @ log_base(nzrvals) - rvals @ P @ log_base(svals)
+    # the relative entropy is guaranteed to be >= 0, so we clamp the
+    # calculated value to 0 to avoid small violations of the lower bound.
+    return max(0, S)
 
 
 def entropy_conditional(rho, selB, base=e, sparse=False):

qutip/mesolve.py

@@ -293,7 +293,17 @@ def _mesolve_QobjEvo(H, c_ops, tlist, args, opt):
     else:
         L_td = H_td
     for op in c_ops:
-        op_td = QobjEvo(op, args, tlist=tlist)
+        # We want to avoid passing tlist where it isn't necessary, to allow a
+        # Hamiltonian/Liouvillian which already _has_ time-dependence not equal
+        # to the mesolve evaluation times to be used in conjunction with
+        # time-independent c_ops.  If we _always_ pass it, it may appear to
+        # QobjEvo that there is a tlist mismatch, even though it is not used.
+        if isinstance(op, Qobj):
+            op_td = QobjEvo(op)
+        elif isinstance(op, QobjEvo):
+            op_td = QobjEvo(op, args)
+        else:
+            op_td = QobjEvo(op, args, tlist=tlist)
         if not issuper(op_td.cte):
             op_td = lindblad_dissipator(op_td)
         L_td += op_td