Clean semidefinite.py

doc/changes/2138.misc

@@ -0,0 +1 @@
+Clean semidefinite

qutip/core/metrics.py

@@ -540,16 +540,10 @@ def dnorm(A, B=None, solver="CVXOPT", verbose=False, force_solve=False,
     # Assume square...
     dim = np.prod(J.dims[0][0])
 
-    # The constraints only depend on the dimension, so
-    # we can cache them efficiently.
-    problem, Jr, Ji, *_ = dnorm_problem(dim)
-
     # Load the parameters with the Choi matrix passed in.
     J_dat = _data.to('csr', J.data).as_scipy()
 
     if not sparse:
-        # The parameters and constraints only depend on the dimension, so
-        # we can cache them efficiently.
         problem, Jr, Ji = dnorm_problem(dim)
 
         # Load the parameters with the Choi matrix passed in.
@@ -561,9 +555,6 @@ def dnorm(A, B=None, solver="CVXOPT", verbose=False, force_solve=False,
                                   J_dat.indptr),
                                  shape=J_dat.shape).toarray()
     else:
-
-        # The parameters do not depend solely on the dimension,
-        # so we can not cache them efficiently.
         problem = dnorm_sparse_problem(dim, J_dat)
 
     problem.solve(solver=solver, verbose=verbose)

qutip/core/operators.py

@@ -6,7 +6,7 @@
 __all__ = ['jmat', 'spin_Jx', 'spin_Jy', 'spin_Jz', 'spin_Jm', 'spin_Jp',
            'spin_J_set', 'sigmap', 'sigmam', 'sigmax', 'sigmay', 'sigmaz',
            'destroy', 'create', 'qeye', 'identity', 'position', 'momentum',
-           'num', 'squeeze', 'squeezing', 'displace', 'commutator',
+           'num', 'squeeze', 'squeezing', 'swap', 'displace', 'commutator',
            'qutrit_ops', 'qdiags', 'phase', 'qzero', 'enr_destroy',
            'enr_identity', 'charge', 'tunneling', 'qft']
 
@@ -1071,3 +1071,31 @@ def qft(dimensions, *, dtype="dense"):
     L, M = np.meshgrid(arr, arr)
     data = np.exp(phase * (L * M)) / np.sqrt(N2)
     return Qobj(data, dims=dimensions).to(dtype)
+
+
+def swap(N, M, *, dtype=None):
+    """
+    Operator that exchanges the order of tensored spaces:
+
+        swap(N, M) @ tensor(ketN, ketM) == tensor(ketM, ketN)
+
+    parameters
+    ----------
+    N : int
+        Number of basis states in the first Hilbert space.
+
+    M : int
+        Number of basis states in the second Hilbert space.
+    """
+    dtype = dtype or settings.core["default_dtype"] or _data.CSR
+
+    if N == 1 and M == 1:
+        return qeye([1, 1], dtype=dtype)
+
+    data = np.ones(N * M)
+    rows = np.arange(N * M + 1)  # last entry is nnz
+    cols = np.ravel(M * np.arange(N)[None, :] + np.arange(M)[:, None])
+    return Qobj(
+        _data.CSR((data, cols, rows), (N * M, N * M)),
+        dims=[[M, N], [N, M]]
+    ).to(dtype)

qutip/core/semidefinite.py

@@ -3,61 +3,47 @@
 """
 This module implements internal-use functions for semidefinite programming.
 """
-
-__all__ = [
-    'dnorm_problem',
-    'dnorm_sparse_problem',
-]
-
 import collections
 import functools
 
 import numpy as np
 import scipy.sparse as sp
+
 # Conditionally import CVXPY
 try:
     import cvxpy
+
+    __all__ = ["dnorm_problem", "dnorm_sparse_problem"]
 except ImportError:
     cvxpy = None
+    __all__ = []
 
-from .tensor import tensor_swap
-from .operators import qeye
-from ..logging_utils import get_logger
-logger = get_logger('qutip.core.semidefinite')
+from .operators import swap
 
-Complex = collections.namedtuple('Complex', ['re', 'im'])
+Complex = collections.namedtuple("Complex", ["re", "im"])
 
 
-def complex_var(rows=1, cols=1, name=None):
+def _complex_var(rows=1, cols=1, name=None):
     return Complex(
         re=cvxpy.Variable((rows, cols), name=(name + "_re") if name else None),
         im=cvxpy.Variable((rows, cols), name=(name + "_im") if name else None),
     )
 
 
-def herm(*Xs):
-    return sum([[X.re == X.re.T, X.im == -X.im.T] for X in Xs], [])
-
-
-def pos_noherm(*Xs):
-    return [
-        cvxpy.bmat([
-            [X.re, -X.im],
-            [X.im, X.re]
-        ]) >> 0
-        for X in Xs
-    ]
-
-
-def pos(*Xs):
-    return pos_noherm(*Xs) + herm(*Xs)
-
-
-def dens(*rhos):
-    return pos(*rhos) + [
-        cvxpy.trace(rho.re) == 1
-        for rho in rhos
+def _make_constraints(*rhos):
+    """
+    Create constraints to ensure definied density operators.
+    """
+    # rhos traces are 1
+    constraints = [cvxpy.trace(rho.re) == 1 for rho in rhos]
+    # rhos are Hermitian
+    for rho in rhos:
+        constraints += [rho.re == rho.re.T] + [rho.im == -rho.im.T]
+    # Non negative
+    constraints += [
+        cvxpy.bmat([[rho.re, -rho.im], [rho.im, rho.re]]) >> 0 for rho in rhos
     ]
+    return constraints
 
 
 def _arr_to_complex(A):
@@ -66,7 +52,7 @@ def _arr_to_complex(A):
     return Complex(re=A, im=np.zeros_like(A))
 
 
-def kron(A, B):
+def _kron(A, B):
     if isinstance(A, np.ndarray):
         A = _arr_to_complex(A)
     if isinstance(B, np.ndarray):
@@ -78,98 +64,68 @@ def kron(A, B):
     )
 
 
-def conj(W, A):
+def _conj(W, A):
     U, V = W.re, W.im
     A, B = A.re, A.im
     return Complex(
         re=(U @ A @ U.T - U @ B @ V.T - V @ A @ V.T - V @ B @ U.T),
-        im=(U @ A @ V.T + U @ B @ U.T + V @ A @ U.T - V @ B @ V.T)
+        im=(U @ A @ V.T + U @ B @ U.T + V @ A @ U.T - V @ B @ V.T),
     )
 
 
-def bmat(B):
-    return Complex(
-        re=cvxpy.bmat([[element.re for element in row] for row in B]),
-        im=cvxpy.bmat([[element.re for element in row] for row in B]),
-    )
-
-
-def dag(X):
-    return Complex(re=X.re.T, im=-X.im.T)
-
-
-def memoize(fn):
-    cache = {}
-
-    @functools.wraps(fn)
-    def memoized(*args):
-        if args in cache:
-            return cache[args]
-        else:
-            ret = fn(*args)
-            cache[args] = ret
-            return ret
-
-    memoized.reset_cache = cache.clear
-    return memoized
-
-
-def qudit_swap(dim):
-    # We should likely generalize this and include it in qip.gates.
-    W = qeye([dim, dim])
-    return tensor_swap(W, (0, 1))
-
-
-@memoize
+@functools.lru_cache
 def initialize_constraints_on_dnorm_problem(dim):
     # Start assembling constraints and variables.
     constraints = []
 
     # Make a complex variable for X.
-    X = complex_var(dim ** 2, dim ** 2, "X")
+    X = _complex_var(dim**2, dim**2, "X")
 
     # Make complex variables for rho0 and rho1.
-    rho0 = complex_var(dim, dim, "rho0")
-    rho1 = complex_var(dim, dim, "rho1")
-    constraints += dens(rho0, rho1)
+    rho0 = _complex_var(dim, dim, "rho0")
+    rho1 = _complex_var(dim, dim, "rho1")
+    constraints += _make_constraints(rho0, rho1)
 
     # Finally, add the tricky positive semidefinite constraint.
     # Since we're using column-stacking, but Watrous used row-stacking,
     # we need to swap the order in Rho0 and Rho1. This is not straightforward,
     # as CVXPY requires that the constant be the first argument. To solve this,
     # We conjugate by SWAP.
-    W = qudit_swap(dim).full()
+    W = swap(dim, dim).full()
     W = Complex(re=W.real, im=W.imag)
-    Rho0 = conj(W, kron(np.eye(dim), rho0))
-    Rho1 = conj(W, kron(np.eye(dim), rho1))
-
-    Y = cvxpy.bmat([
-        [Rho0.re, X.re,      -Rho0.im, -X.im],
-        [X.re.T, Rho1.re,    X.im.T, -Rho1.im],
-
-        [Rho0.im, X.im,      Rho0.re, X.re],
-        [-X.im.T, Rho1.im,   X.re.T, Rho1.re],
-    ])
+    Rho0 = _conj(W, _kron(np.eye(dim), rho0))
+    Rho1 = _conj(W, _kron(np.eye(dim), rho1))
+
+    Y = cvxpy.bmat(
+        [
+            [Rho0.re, X.re, -Rho0.im, -X.im],
+            [X.re.T, Rho1.re, X.im.T, -Rho1.im],
+            [Rho0.im, X.im, Rho0.re, X.re],
+            [-X.im.T, Rho1.im, X.re.T, Rho1.re],
+        ]
+    )
     constraints += [Y >> 0]
 
-    logger.debug("Using %d constraints.", len(constraints))
-
     return X, constraints
 
 
 def dnorm_problem(dim):
+    """
+    Creade the cvxpy ``Problem`` for the dnorm metric using dense arrays
+    """
     X, constraints = initialize_constraints_on_dnorm_problem(dim)
     Jr = cvxpy.Parameter((dim**2, dim**2))
     Ji = cvxpy.Parameter((dim**2, dim**2))
     # The objective, however, depends on J.
-    objective = cvxpy.Maximize(cvxpy.trace(
-        Jr.T @ X.re + Ji.T @ X.im
-    ))
+    objective = cvxpy.Maximize(cvxpy.trace(Jr.T @ X.re + Ji.T @ X.im))
     problem = cvxpy.Problem(objective, constraints)
     return problem, Jr, Ji
 
 
 def dnorm_sparse_problem(dim, J_dat):
+    """
+    Creade the cvxpy ``Problem`` for the dnorm metric using sparse arrays
+    """
     X, constraints = initialize_constraints_on_dnorm_problem(dim)
     J_val = J_dat.tocoo()
 
@@ -186,10 +142,11 @@ def adapt_sparse_params(A_val, dim):
         A_cols = np.arange(A_nnz.size)
         # We are pushing the data on the location of the nonzero elements
         # to the nonzero rows of A_indexer
-        A_Indexer = sp.coo_matrix((A_data, (A_rows, A_cols)),
-                                  shape=(side_size**2, A_nnz.size))
+        A_Indexer = sp.coo_matrix(
+            (A_data, (A_rows, A_cols)), shape=(side_size**2, A_nnz.size)
+        )
         # We get finaly the sparse matrix A which we wanted
-        A = cvxpy.reshape(A_Indexer @ A_nnz, (side_size, side_size), order='C')
+        A = cvxpy.reshape(A_Indexer @ A_nnz, (side_size, side_size), order="C")
         A_nnz.value = A_val.data
         return A
 
@@ -200,9 +157,7 @@ def adapt_sparse_params(A_val, dim):
     Ji = adapt_sparse_params(Ji_val, dim)
 
     # The objective, however, depends on J.
-    objective = cvxpy.Maximize(cvxpy.trace(
-        Jr.T @ X.re + Ji.T @ X.im
-    ))
+    objective = cvxpy.Maximize(cvxpy.trace(Jr.T @ X.re + Ji.T @ X.im))
 
     problem = cvxpy.Problem(objective, constraints)
     return problem

qutip/tests/core/test_operators.py

@@ -293,3 +293,12 @@ def test_qft(dims):
         fft = np.fft.fft(qft[:,i])
         fft /= np.sum(fft)
         np.testing.assert_allclose(fft, target, atol=1e-16 * N)
+
+
+@pytest.mark.parametrize('N', [1, 3, 5, 8])
+@pytest.mark.parametrize('M', [1, 3, 5, 8])
+def test_swap(N, M):
+    ket1 = qutip.rand_ket(N)
+    ket2 = qutip.rand_ket(M)
+
+    assert qutip.swap(N, M) @ (ket1 & ket2) == (ket2 & ket1)