diff --git a/.github/workflows/tests.yml b/.github/workflows/tests.yml index 6f89efa2a4..eceee8cc60 100644 --- a/.github/workflows/tests.yml +++ b/.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 diff --git a/qutip.bib b/CITATION.bib similarity index 100% rename from qutip.bib rename to CITATION.bib diff --git a/doc/apidoc/functions.rst b/doc/apidoc/functions.rst index d13de7f3a6..801e89b010 100644 --- a/doc/apidoc/functions.rst +++ b/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: diff --git a/doc/conf.py b/doc/conf.py index 00b846a2b5..dd6fb9c896 100755 --- a/doc/conf.py +++ b/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. diff --git a/doc/contributors.rst b/doc/contributors.rst index 3a6e2f3494..1a2eba701a 100644 --- a/doc/contributors.rst +++ b/doc/contributors.rst @@ -79,7 +79,7 @@ Contributors - Louis Tessler - Lucas Verney - Marco David -- Marek +- Marek Narozniak - Markus Baden - Martín Sande - Mateo Laguna diff --git a/doc/guide/dynamics/dynamics-floquet.rst b/doc/guide/dynamics/dynamics-floquet.rst index cc9b7420e0..95976817d4 100644 --- a/doc/guide/dynamics/dynamics-floquet.rst +++ b/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: diff --git a/doc/guide/dynamics/dynamics-monte.rst b/doc/guide/dynamics/dynamics-monte.rst index d98e420797..db3e4779ed 100644 --- a/doc/guide/dynamics/dynamics-monte.rst +++ b/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. diff --git a/doc/requirements.txt b/doc/requirements.txt index f6da069ab6..58c3f66db5 100644 --- a/doc/requirements.txt +++ b/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 diff --git a/doc/static/site.css b/doc/static/site.css index 8c6988ea63..2a0a81fca1 100644 --- a/doc/static/site.css +++ b/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; } diff --git a/qutip/control/dynamics.py b/qutip/control/dynamics.py index 0efe66e606..d5602aef16 100644 --- a/qutip/control/dynamics.py +++ b/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,7 +1616,7 @@ 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 @@ -1638,6 +1624,8 @@ 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].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): diff --git a/qutip/control/propcomp.py b/qutip/control/propcomp.py index 90a9e305ea..a7e4c8d9f9 100644 --- a/qutip/control/propcomp.py +++ b/qutip/control/propcomp.py @@ -306,16 +306,11 @@ def _get_aug_mat(self, k, j): E = dyn._get_phased_ctrl_dyn_gen(k, j).data*dyn.tau[k] Z = sp.csr_matrix(dg.data.shape) aug = Qobj(sp.vstack([sp.hstack([A, E]), sp.hstack([Z, A])])) - elif dyn.oper_dtype == np.ndarray: + else: A = dg*dyn.tau[k] E = dyn._get_phased_ctrl_dyn_gen(k, j)*dyn.tau[k] Z = np.zeros(dg.shape) aug = np.vstack([np.hstack([A, E]), np.hstack([Z, A])]) - else: - A = dg*dyn.tau[k] - E = dyn._get_phased_ctrl_dyn_gen(k, j)*dyn.tau[k] - Z = dg*0.0 - aug = sp.vstack([sp.hstack([A, E]), sp.hstack([Z, A])]) return aug def _compute_prop_grad(self, k, j, compute_prop=True): @@ -388,7 +383,7 @@ def _compute_prop_grad(self, k, j, compute_prop=True): prop_grad_dense = la.expm_frechet(A, E, compute_expm=False) prop_grad = Qobj(prop_grad_dense, dims=dyn.dyn_dims) - elif dyn.oper_dtype == np.ndarray: + else: A = dyn._get_phased_dyn_gen(k)*dyn.tau[k] E = dyn._get_phased_ctrl_dyn_gen(k, j)*dyn.tau[k] if compute_prop: @@ -396,19 +391,6 @@ def _compute_prop_grad(self, k, j, compute_prop=True): else: prop_grad = la.expm_frechet(A, E, compute_expm=False) - else: - # Assuming some sparse matrix - spcls = dyn._dyn_gen[k].__class__ - A = (dyn._get_phased_dyn_gen(k)*dyn.tau[k]).toarray() - E = (dyn._get_phased_ctrl_dyn_gen(k, j)*dyn.tau[k]).toarray() - if compute_prop: - prop_dense, prop_grad_dense = la.expm_frechet(A, E) - prop = spcls(prop_dense) - prop_grad = spcls(prop_grad_dense) - else: - prop_grad_dense = la.expm_frechet(A, E, compute_expm=False) - prop_grad = spcls(prop_grad_dense) - if compute_prop: return prop, prop_grad else: diff --git a/qutip/cy/mcsolve.pyx b/qutip/cy/mcsolve.pyx index 2e3e1a5db5..5b11f81d98 100644 --- a/qutip/cy/mcsolve.pyx +++ b/qutip/cy/mcsolve.pyx @@ -32,15 +32,15 @@ cdef np.ndarray[complex, ndim=1] normalize(complex[::1] psi): cdef class CyMcOde: cdef: - int steady_state, store_states, col_args + int steady_state, store_states, col_args, store_final_state int norm_steps, l_vec, num_ops double norm_t_tol, norm_tol list collapses list collapses_args list c_ops list n_ops - complex[:,::1] states_out - complex[:,::1] ss_out + object states_out + object ss_out double[::1] n_dp def __init__(self, ss, opt): @@ -51,6 +51,7 @@ cdef class CyMcOde: self.norm_tol = opt.norm_tol self.steady_state = opt.steady_state_average self.store_states = opt.store_states or opt.average_states + self.store_final_state = opt.store_final_state self.collapses = [] self.l_vec = self.c_ops[0].cte.shape[0] self.num_ops = len(ss.td_n_ops) @@ -77,10 +78,11 @@ cdef class CyMcOde: @cython.wraparound(False) cdef void sumsteadystate(self, complex[::1] state): cdef int ii, jj, l_vec + cdef complex [:, ::1] _ss_out = self.ss_out l_vec = state.shape[0] for ii in range(l_vec): - for jj in range(l_vec): - self.ss_out[ii,jj] += state[ii]*conj(state[jj]) + for jj in range(l_vec): + _ss_out[ii,jj] += state[ii] * conj(state[jj]) @cython.boundscheck(False) @@ -94,6 +96,7 @@ cdef class CyMcOde: cdef int ii, which, k cdef double norm2_prev, norm2_psi cdef double t_prev + cdef complex [:, ::1] _states_out if self.steady_state: self.sumsteadystate(out_psi) @@ -102,8 +105,11 @@ cdef class CyMcOde: self.states_out = np.zeros((num_times, self.l_vec), dtype=complex) for ii in range(self.l_vec): self.states_out[0, ii] = out_psi[ii] - else: + elif self.store_final_state: self.states_out = np.zeros((1, self.l_vec), dtype=complex) + else: + self.states_out = None + _states_out = self.states_out e_call.step(0, out_psi) rand_vals = prng.rand(2) @@ -152,11 +158,11 @@ cdef class CyMcOde: self.sumsteadystate(out_psi) if self.store_states: for ii in range(self.l_vec): - self.states_out[k, ii] = out_psi[ii] - if not self.store_states: + _states_out[k, ii] = out_psi[ii] + if not self.store_states and self.store_final_state: for ii in range(self.l_vec): - self.states_out[0, ii] = out_psi[ii] - return np.array(self.states_out), np.array(self.ss_out), self.collapses + _states_out[0, ii] = out_psi[ii] + return self.states_out, self.ss_out, self.collapses @cython.cdivision(True) @@ -341,8 +347,10 @@ cdef class CyMcOdeDiag(CyMcOde): self.states_out = np.zeros((num_times, self.l_vec), dtype=complex) for ii in range(self.l_vec): self.states_out[0, ii] = out_psi[ii] - else: + elif self.store_final_state: self.states_out = np.zeros((1, self.l_vec), dtype=complex) + else: + self.states_out = None e_call.step(0, out_psi) rand_vals = prng.rand(2) @@ -367,10 +375,10 @@ cdef class CyMcOdeDiag(CyMcOde): if self.store_states: for ii in range(self.l_vec): self.states_out[k, ii] = out_psi[ii] - if not self.store_states: + if not self.store_states and self.store_final_state: for ii in range(self.l_vec): self.states_out[0, ii] = out_psi[ii] - return np.array(self.states_out), np.array(self.ss_out), self.collapses + return self.states_out, self.ss_out, self.collapses @cython.cdivision(True) @cython.boundscheck(False) diff --git a/qutip/mesolve.py b/qutip/mesolve.py index 180ac569a3..84b7389ded 100644 --- a/qutip/mesolve.py +++ b/qutip/mesolve.py @@ -127,9 +127,13 @@ def mesolve(H, rho0, tlist, c_ops=None, e_ops=None, args=None, options=None, single collapse operator, or list of collapse operators, or a list of Liouvillian superoperators. - e_ops : None / list of :class:`qutip.Qobj` / callback function single - single operator or list of operators for which to evaluate - expectation values. + e_ops : None / list / callback function, optional + A list of operators as `Qobj` and/or callable functions (can be mixed) + or a single callable function. For operators, the result's expect will + be computed by :func:`qutip.expect`. For callable functions, they are + called as ``f(t, state)`` and return the expectation value. + A single callback's expectation value can be any type, but a callback + part of a list must return a number as the expectation value. args : None / *dictionary* dictionary of parameters for time-dependent Hamiltonians and @@ -458,6 +462,9 @@ def _generic_ode_solve(func, ode_args, rho0, tlist, e_ops, opt, opt.store_states = True else: for op in e_ops: + if not isinstance(op, Qobj) and callable(op): + output.expect.append(np.zeros(n_tsteps, dtype=complex)) + continue e_ops_data.append(spre(op).data) if op.isherm and rho0.isherm: output.expect.append(np.zeros(n_tsteps)) @@ -504,6 +511,9 @@ def get_curr_state_data(r): output.expect.append(e_ops(t, rho_t)) for m in range(n_expt_op): + if not isinstance(e_ops[m], Qobj) and callable(e_ops[m]): + output.expect[m][t_idx] = e_ops[m](t, rho_t) + continue output.expect[m][t_idx] = expect_rho_vec(e_ops_data[m], r.y, e_ops[m].isherm and rho0.isherm) diff --git a/qutip/operators.py b/qutip/operators.py index 7d787a176c..7d5e011ce6 100644 --- a/qutip/operators.py +++ b/qutip/operators.py @@ -859,7 +859,7 @@ def enr_destroy(dims, excitations): a_ops = [sp.lil_matrix((nstates, nstates), dtype=np.complex128) for _ in range(len(dims))] - for n1, state1 in idx2state.items(): + for n1, state1 in enumerate(idx2state): for idx, s in enumerate(state1): # if s > 0, the annihilation operator of mode idx has a non-zero # entry with one less excitation in mode idx in the final state diff --git a/qutip/qip/circuit.py b/qutip/qip/circuit.py index e92127e713..b3be05a0ea 100644 --- a/qutip/qip/circuit.py +++ b/qutip/qip/circuit.py @@ -364,6 +364,8 @@ class QubitCircuit: A list of integer for the dimension of each composite system. e.g [2,2,2,2,2] for 5 qubits system. If None, qubits system will be the default option. + num_cbits : int + Number of classical bits in the system. Examples -------- diff --git a/qutip/qip/operations/gates.py b/qutip/qip/operations/gates.py index 7c72cf66dd..7b58b9bdf9 100644 --- a/qutip/qip/operations/gates.py +++ b/qutip/qip/operations/gates.py @@ -892,11 +892,9 @@ def hadamard_transform(N=1): Quantum object representation of the N-qubit Hadamard gate. """ - data = 2 ** (-N / 2) * np.array([[(-1) ** _hamming_distance(i & j) - for i in range(2 ** N)] - for j in range(2 ** N)]) + H = Qobj([[1, 1], [1, -1]]) / np.sqrt(2) - return Qobj(data, dims=[[2] * N, [2] * N]) + return tensor([H] * N) def _flatten(lst): diff --git a/qutip/qobj.py b/qutip/qobj.py index 8d280d6b18..bf82a4a665 100644 --- a/qutip/qobj.py +++ b/qutip/qobj.py @@ -693,7 +693,7 @@ def __pow__(self, n, m=None): # calculates powers of Qobj """ POWER operation. """ - if self.type not in ['oper', 'super']: + if self.shape[0] != self.shape[1]: raise Exception("Raising a qobj to some power works only for " + "operators and super-operators (square matrices).") diff --git a/qutip/sesolve.py b/qutip/sesolve.py index 75522d36c5..2374f2cc91 100644 --- a/qutip/sesolve.py +++ b/qutip/sesolve.py @@ -54,9 +54,14 @@ def sesolve(H, psi0, tlist, e_ops=None, args=None, options=None, tlist : array_like of float List of times for :math:`t`. - e_ops : list of :class:`~Qobj` or callback function, optional - Single operator or list of operators for which to evaluate expectation - values. For operator evolution, the overlap is computed: :: + e_ops : None / list / callback function, optional + A list of operators as `Qobj` and/or callable functions (can be mixed) + or a single callable function. For callable functions, they are called + as ``f(t, state)`` and return the expectation value. A single + callback's expectation value can be any type, but a callback part of a + list must return a number as the expectation value. For operators, the + result's expect will be computed by :func:`qutip.expect` when the state + is a ``ket``. For operator evolution, the overlap is computed by: :: (e_ops[i].dag() * op(t)).tr() @@ -275,16 +280,25 @@ def _generic_ode_solve(func, ode_args, psi0, tlist, e_ops, opt, opt.store_states = True else: for op in e_ops: + if not isinstance(op, Qobj) and callable(op): + output.expect.append(np.zeros(n_tsteps, dtype=complex)) + continue if op.isherm: output.expect.append(np.zeros(n_tsteps)) else: output.expect.append(np.zeros(n_tsteps, dtype=complex)) if oper_evo: for e in e_ops: - e_ops_data.append(e.dag().data) + if isinstance(e, Qobj): + e_ops_data.append(e.dag().data) + continue + e_ops_data.append(e) else: for e in e_ops: - e_ops_data.append(e.data) + if isinstance(e, Qobj): + e_ops_data.append(e.data) + continue + e_ops_data.append(e) else: raise TypeError("Expectation parameter must be a list or a function") @@ -311,7 +325,8 @@ def get_curr_state_data(r): "the allowed number of substeps by increasing " "the nsteps parameter in the Options class.") # get the current state / oper data if needed - if opt.store_states or opt.normalize_output or n_expt_op > 0 or expt_callback: + if opt.store_states or opt.normalize_output \ + or n_expt_op > 0 or expt_callback: cdata = get_curr_state_data(r) if opt.normalize_output: @@ -343,9 +358,17 @@ def get_curr_state_data(r): if oper_evo: for m in range(n_expt_op): + if callable(e_ops_data[m]): + func = e_ops_data[m] + output.expect[m][t_idx] = func(t, Qobj(cdata, dims=dims)) + continue output.expect[m][t_idx] = (e_ops_data[m] * cdata).trace() else: for m in range(n_expt_op): + if callable(e_ops_data[m]): + func = e_ops_data[m] + output.expect[m][t_idx] = func(t, Qobj(cdata, dims=dims)) + continue output.expect[m][t_idx] = cy_expect_psi(e_ops_data[m], cdata, e_ops[m].isherm) diff --git a/qutip/simdiag.py b/qutip/simdiag.py index 505738cc27..6b57f82cf4 100644 --- a/qutip/simdiag.py +++ b/qutip/simdiag.py @@ -5,7 +5,43 @@ from qutip.qobj import Qobj -def simdiag(ops, evals=True): +def _degen(tol, vecs, ops, i=0): + """ + Private function that finds eigen vals and vecs for degenerate matrices.. + """ + if len(ops) == i: + return vecs + + # New eigenvectors are sometime not orthogonal. + for j in range(1, vecs.shape[1]): + for k in range(j): + dot = vecs[:, j].dot(vecs[:, k].conj()) + if np.abs(dot) > tol: + vecs[:, j] = ((vecs[:, j] - dot * vecs[:, k]) + / (1 - np.abs(dot)**2)**0.5) + + subspace = vecs.conj().T @ ops[i].data @ vecs + eigvals, eigvecs = la.eig(subspace) + + perm = np.argsort(eigvals) + eigvals = eigvals[perm] + + vecs_new = vecs @ eigvecs[:, perm] + for k in range(len(eigvals)): + vecs_new[:, k] = vecs_new[:, k] / la.norm(vecs_new[:, k]) + + k = 0 + while k < len(eigvals): + ttol = max(tol, tol * abs(eigvals[k])) + inds, = np.where(abs(eigvals - eigvals[k]) < ttol) + if len(inds) > 1: # if at least 2 eigvals are degenerate + vecs_new[:, inds] = _degen(tol, vecs_new[:, inds], ops, i+1) + k = inds[-1] + 1 + return vecs_new + + +def simdiag(ops, evals: bool = True, *, + tol: float = 1e-14, safe_mode: bool = True): """Simultaneous diagonalization of commuting Hermitian matrices. Parameters @@ -14,6 +50,18 @@ def simdiag(ops, evals=True): ``list`` or ``array`` of qobjs representing commuting Hermitian operators. + evals : bool [True] + Whether to return the eigenvalues for each ops and eigenvectors or just + the eigenvectors. + + tol : float [1e-14] + Tolerance for detecting degenerate eigenstates. + + safe_mode : bool [True] + Whether to check that all ops are Hermitian and commuting. If set to + ``False`` and operators are not commuting, the eigenvectors returned + will often be eigenvectors of only the first operator. + Returns -------- eigs : tuple @@ -22,21 +70,16 @@ def simdiag(ops, evals=True): operator. """ - tol = 1e-14 - start_flag = 0 - if not any(ops): - raise ValueError('Need at least one input operator.') - if not isinstance(ops, (list, np.ndarray)): - ops = np.array([ops]) - num_ops = len(ops) + if not ops: + raise ValueError("No input matrices.") + N = ops[0].shape[0] + num_ops = len(ops) if safe_mode else 0 for jj in range(num_ops): A = ops[jj] shape = A.shape if shape[0] != shape[1]: raise TypeError('Matricies must be square.') - if start_flag == 0: - s = shape[0] - if s != shape[0]: + if shape[0] != N: raise TypeError('All matrices. must be the same shape') if not A.isherm: raise TypeError('Matricies must be Hermitian') @@ -45,75 +88,34 @@ def simdiag(ops, evals=True): if (A * B - B * A).norm() / (A * B).norm() > tol: raise TypeError('Matricies must commute.') - A = ops[0] - eigvals, eigvecs = la.eig(A.full()) - zipped = list(zip(-eigvals, range(len(eigvals)))) - zipped.sort() - ds, perm = zip(*zipped) - ds = -np.real(np.array(ds)) - perm = np.array(perm) - eigvecs_array = np.array( - [np.zeros((A.shape[0], 1), dtype=complex) for k in range(A.shape[0])]) - - for kk in range(len(perm)): # matrix with sorted eigvecs in columns - eigvecs_array[kk][:, 0] = eigvecs[:, perm[kk]] + eigvals, eigvecs = la.eigh(ops[0].full()) + perm = np.argsort(eigvals) + eigvecs = eigvecs[:, perm] + eigvals = eigvals[perm] + k = 0 - rng = np.arange(len(eigvals)) - while k < len(ds): + while k < N: # find degenerate eigenvalues, get indicies of degenerate eigvals - inds = np.array(abs(ds - ds[k]) < max(tol, tol * abs(ds[k]))) - inds = rng[inds] + ttol = max(tol, tol * abs(eigvals[k])) + inds, = np.where(abs(eigvals - eigvals[k]) < ttol) if len(inds) > 1: # if at least 2 eigvals are degenerate - eigvecs_array[inds] = degen(tol, eigvecs_array[inds], ops[1:]) - k = max(inds) + 1 - eigvals_out = np.zeros((num_ops, len(ds)), dtype=np.float64) - kets_out = np.empty((len(ds),), dtype=object) - kets_out[:] = [ - Qobj(eigvecs_array[j] / la.norm(eigvecs_array[j]), + eigvecs[:, inds] = _degen(tol, eigvecs[:, inds], ops, 1) + k = inds[-1] + 1 + + for k in range(N): + eigvecs[:, k] = eigvecs[:, k] / la.norm(eigvecs[:, k]) + + kets_out = [ + Qobj(eigvecs[:, j], dims=[ops[0].dims[0], [1]], shape=[ops[0].shape[0], 1]) - for j in range(len(ds)) + for j in range(N) ] + eigvals_out = np.zeros((len(ops), N), dtype=np.float64) if not evals: return kets_out else: - for kk in range(num_ops): - for j in range(len(ds)): - eigvals_out[kk, j] = np.real(np.dot( - eigvecs_array[j].conj().T, - ops[kk].data * eigvecs_array[j])) + for kk in range(len(ops)): + for j in range(N): + eigvals_out[kk, j] = ops[kk].matrix_element(kets_out[j], + kets_out[j]).real return eigvals_out, kets_out - - -def degen(tol, in_vecs, ops): - """ - Private function that finds eigen vals and vecs for degenerate matrices.. - """ - n = len(ops) - if n == 0: - return in_vecs - A = ops[0] - vecs = np.column_stack(in_vecs) - eigvals, eigvecs = la.eig(np.dot(vecs.conj().T, A.data.dot(vecs))) - zipped = list(zip(-eigvals, range(len(eigvals)))) - zipped.sort() - ds, perm = zip(*zipped) - ds = -np.real(np.array(ds)) - perm = np.array(perm) - vecsperm = np.zeros(eigvecs.shape, dtype=complex) - for kk in range(len(perm)): # matrix with sorted eigvecs in columns - vecsperm[:, kk] = eigvecs[:, perm[kk]] - vecs_new = np.dot(vecs, vecsperm) - vecs_out = np.array( - [np.zeros((A.shape[0], 1), dtype=complex) for k in range(len(ds))]) - for kk in range(len(perm)): # matrix with sorted eigvecs in columns - vecs_out[kk][:, 0] = vecs_new[:, kk] - k = 0 - rng = np.arange(len(ds)) - while k < len(ds): - inds = np.array(abs(ds - ds[k]) < max( - tol, tol * abs(ds[k]))) # get indicies of degenerate eigvals - inds = rng[inds] - if len(inds) > 1: # if at least 2 eigvals are degenerate - vecs_out[inds] = degen(tol, vecs_out[inds], ops[1:n]) - k = max(inds) + 1 - return vecs_out diff --git a/qutip/solver.py b/qutip/solver.py index 67634d4b7b..307db3911e 100644 --- a/qutip/solver.py +++ b/qutip/solver.py @@ -43,14 +43,26 @@ def __init__(self, e_ops=[], super_=False): self.e_ops = e_ops if isinstance(e_ops, list): self.e_num = len(e_ops) - self.e_ops_isherm = [e.isherm for e in e_ops] - if not super_: - self.e_ops_qoevo = np.array([QobjEvo(e) for e in e_ops], - dtype=object) - else: - self.e_ops_qoevo = np.array([QobjEvo(spre(e)) for e in e_ops], - dtype=object) - [op.compile() for op in self.e_ops_qoevo] + e_ops_qoevo = [] + e_ops_isherm = [] + for e in e_ops: + if isinstance(e, (Qobj, QobjEvo)): + e_ops_isherm.append(e.isherm) + e_ops_qoevo_entry = None + if not super_: + e_ops_qoevo_entry = QobjEvo(e) + else: + e_ops_qoevo_entry = QobjEvo(spre(e)) + e_ops_qoevo_entry.compile() + e_ops_qoevo.append(e_ops_qoevo_entry) + elif callable(e): + e_ops_isherm.append(None) + e_ops_qoevo.append(e) + else: + raise TypeError("Expectation value list entry needs to be " + "either a function either an operator") + self.e_ops_isherm = e_ops_isherm + self.e_ops_qoevo = np.array(e_ops_qoevo, dtype=object) elif callable(e_ops): self.isfunc = True self.e_num = 1 @@ -79,8 +91,12 @@ def step(self, iter_, state): else: t = self.tlist[iter_] for ii in range(self.e_num): - self.raw_out[ii, iter_] = \ - self.e_ops_qoevo[ii].compiled_qobjevo.expect(t, state) + if isinstance(self.e_ops_qoevo[ii], QobjEvo): + self.raw_out[ii, iter_] = \ + self.e_ops_qoevo[ii].compiled_qobjevo.expect(t, state) + elif callable(self.e_ops_qoevo[ii]): + self.raw_out[ii, iter_] = \ + self.e_ops_qoevo[ii](t, state) def finish(self): if self.isfunc: diff --git a/qutip/states.py b/qutip/states.py index 5cada3b1f1..10ac9e3296 100644 --- a/qutip/states.py +++ b/qutip/states.py @@ -11,6 +11,7 @@ import numpy as np from numpy import arange, conj, prod import scipy.sparse as sp +import itertools from qutip.qobj import Qobj from qutip.operators import destroy, jmat @@ -715,67 +716,62 @@ def bra(seq, dim=2): # # quantum state number helper functions # -def state_number_enumerate(dims, excitations=None, state=None, idx=0, nexc=0): +def state_number_enumerate(dims, excitations=None): """ - An iterator that enumerate all the state number arrays (quantum numbers on - the form [n1, n2, n3, ...]) for a system with dimensions given by dims. + An iterator that enumerates all the state number tuples (quantum numbers of + the form (n1, n2, n3, ...)) for a system with dimensions given by dims. Example: >>> for state in state_number_enumerate([2,2]): # doctest: +SKIP >>> print(state) # doctest: +SKIP - [ 0 0 ] - [ 0 1 ] - [ 1 0 ] - [ 1 1 ] + ( 0 0 ) + ( 0 1 ) + ( 1 0 ) + ( 1 1 ) Parameters ---------- dims : list or array The quantum state dimensions array, as it would appear in a Qobj. - state : list - Current state in the iteration. Used internally. - excitations : integer (None) Restrict state space to states with excitation numbers below or equal to this value. - idx : integer - Current index in the iteration. Used internally. - - nexc : integer - Number of excitations in modes [0..idx-1]. Used internally. - Returns ------- - state_number : list - Successive state number arrays that can be used in loops and other + state_number : tuple + Successive state number tuples that can be used in loops and other iterations, using standard state enumeration *by definition*. """ - if state is None: - state = np.zeros(len(dims), dtype=int) - - if idx == len(dims): - if excitations is None: - yield np.array(state) - else: - yield tuple(state) - else: - if excitations is None: - nlim = dims[idx] - else: - # modes [0..idx-1] have nexc excitations, - # so mode idx can have at most excitations-nexc excitations - nlim = min(dims[idx], 1 + excitations - nexc) - - for n in range(nlim): - state[idx] = n - for s in state_number_enumerate(dims, excitations, - state, idx + 1, nexc + n): - yield s + if excitations is None: + # in this case, state numbers are a direct product + yield from itertools.product(*(range(d) for d in dims)) + return + + # From here on, excitations is not None + + # General idea of algorithm: add excitations one by one in last mode (idx = + # len(dims)-1), and carry over to the next index when the limit is reached. + # Keep track of the number of excitations while doing so to avoid having to + # do explicit sums over the states. + state = (0,)*len(dims) + nexc = 0 + while True: + yield state + idx = len(dims) - 1 + state = state[:idx] + (state[idx]+1,) + nexc += 1 + while nexc > excitations or state[idx] >= dims[idx]: + # remove all excitations in mode idx, add one in idx-1 + idx -= 1 + if idx < 0: + return + nexc -= state[idx+1] - 1 + state = state[:idx] + (state[idx]+1, 0) + state[idx+2:] def state_number_index(dims, state): @@ -803,8 +799,7 @@ def state_number_index(dims, state): ordering. """ - return int( - sum([state[i] * prod(dims[i + 1:]) for i, d in enumerate(dims)])) + return np.ravel_multi_index(state, dims) def state_index_number(dims, index): @@ -827,20 +822,12 @@ def state_index_number(dims, index): Returns ------- - state : list - The state number array corresponding to index `index` in standard + state : tuple + The state number tuple corresponding to index `index` in standard enumeration ordering. """ - state = np.empty_like(dims) - - D = np.concatenate([np.flipud(np.cumprod(np.flipud(dims[1:]))), [1]]) - - for n in range(len(dims)): - state[n] = index / D[n] - index -= state[n] * D[n] - - return list(state) + return np.unravel_index(index, dims) def state_number_qobj(dims, state): @@ -878,7 +865,8 @@ def state_number_qobj(dims, state): """ - return tensor([fock(dims[i], s) for i, s in enumerate(state)]) + assert len(state) == len(dims) + return tensor([fock(d, s) for d, s in zip(dims, state)]) # @@ -900,19 +888,16 @@ def enr_state_dictionaries(dims, excitations): Returns ------- - nstates, state2idx, idx2state: integer, dict, dict + nstates, state2idx, idx2state: integer, dict, list The number of states `nstates`, a dictionary for looking up state - indices from a state tuple, and a dictionary for looking up state - state tuples from state indices. + indices from a state tuple, and a list containing the state tuples + ordered by state indices. state2idx and idx2state are reverses of + each other, i.e., state2idx[idx2state[idx]] = idx and + idx2state[state2idx[state]] = state. """ - nstates = 0 - state2idx = {} - idx2state = {} - - for state in state_number_enumerate(dims, excitations): - state2idx[state] = nstates - idx2state[nstates] = state - nstates += 1 + idx2state = list(state_number_enumerate(dims, excitations)) + state2idx = {state: idx for idx, state in enumerate(idx2state)} + nstates = len(idx2state) return nstates, state2idx, idx2state @@ -996,8 +981,7 @@ def enr_thermal_dm(dims, excitations, n): else: n = np.asarray(n) - diags = [np.prod((n / (n + 1)) ** np.array(state)) - for idx, state in idx2state.items()] + diags = [np.prod((n / (n + 1)) ** np.array(state)) for state in idx2state] diags /= np.sum(diags) data = sp.spdiags(diags, 0, nstates, nstates, format='csr') diff --git a/qutip/steadystate.py b/qutip/steadystate.py index 3b055c21d9..1f60f942b2 100644 --- a/qutip/steadystate.py +++ b/qutip/steadystate.py @@ -4,8 +4,8 @@ collapse operators. """ -__all__ = ['steadystate', 'steady', 'build_preconditioner', - 'pseudo_inverse'] +__all__ = ['steadystate', 'steady', 'steadystate_floquet', + 'build_preconditioner', 'pseudo_inverse'] import functools import time @@ -242,10 +242,10 @@ def steadystate(A, c_op_list=[], method='direct', solver=None, **kwargs): solver = 'mkl' elif solver == 'mkl' and \ (method not in ['direct', 'power']): - raise Exception('MKL solver only for direct or power methods.') + raise ValueError('MKL solver only for direct or power methods.') elif solver not in ['scipy', 'mkl']: - raise Exception('Invalid solver kwarg.') + raise ValueError('Invalid solver kwarg.') ss_args = _default_steadystate_args() ss_args['method'] = method @@ -258,7 +258,7 @@ def steadystate(A, c_op_list=[], method='direct', solver=None, **kwargs): if key in ss_args.keys(): ss_args[key] = kwargs[key] else: - raise Exception( + raise TypeError( "Invalid keyword argument '"+key+"' passed to steadystate.") # Set column perm to NATURAL if using RCM and not specified by user @@ -914,6 +914,100 @@ def _iter_count(r): return rhoss +def steadystate_floquet(H_0, c_ops, Op_t, w_d=1.0, n_it=3, sparse=False): + """ + Calculates the effective steady state for a driven + system with a time-dependent cosinusoidal term: + + .. math:: + + \\mathcal{\\hat{H}}(t) = \\hat{H}_0 + + \\mathcal{\\hat{O}} \\cos(\\omega_d t) + + Parameters + ---------- + H_0 : :obj:`~Qobj` + A Hamiltonian or Liouvillian operator. + + c_ops : list + A list of collapse operators. + + Op_t : :obj:`~Qobj` + The the interaction operator which is multiplied by the cosine + + w_d : float, default 1.0 + The frequency of the drive + + n_it : int, default 3 + The number of iterations for the solver + + sparse : bool, default False + Solve for the steady state using sparse algorithms. + Actually, dense seems to be faster. + + Returns + ------- + dm : qobj + Steady state density matrix. + + .. note:: + See: Sze Meng Tan, + https://copilot.caltech.edu/documents/16743/qousersguide.pdf, + Section (10.16) + """ + if sparse: + N = H_0.shape[0] + + L_0 = liouvillian(H_0, c_ops).data.tocsc() + L_t = liouvillian(Op_t) + L_p = (0.5 * L_t).data.tocsc() + # L_p and L_m correspond to the positive and negative + # frequency terms respectively. + # They are independent in the model, so we keep both names. + L_m = L_p + L_p_array = L_p.todense() + L_m_array = L_p_array + + Id = sp.eye(N ** 2, format="csc", dtype=np.complex128) + S = T = sp.csc_matrix((N ** 2, N ** 2), dtype=np.complex128) + + for n_i in np.arange(n_it, 0, -1): + L = sp.csc_matrix(L_0 - 1j * n_i * w_d * Id + L_m.dot(S)) + L.sort_indices() + LU = splu(L) + S = - LU.solve(L_p_array) + + L = sp.csc_matrix(L_0 + 1j * n_i * w_d * Id + L_p.dot(T)) + L.sort_indices() + LU = splu(L) + T = - LU.solve(L_m_array) + + M_subs = L_0 + L_m.dot(S) + L_p.dot(T) + else: + N = H_0.shape[0] + + L_0 = liouvillian(H_0, c_ops).full() + L_t = liouvillian(Op_t) + L_p = (0.5 * L_t).full() + L_m = L_p + + Id = np.eye(N ** 2) + S, T = np.zeros((N ** 2, N ** 2)), np.zeros((N ** 2, N ** 2)) + + for n_i in np.arange(n_it, 0, -1): + L = L_0 - 1j * n_i * w_d * Id + np.matmul(L_m, S) + lu, piv = la.lu_factor(L) + S = - la.lu_solve((lu, piv), L_p) + + L = L_0 + 1j * n_i * w_d * Id + np.matmul(L_p, T) + lu, piv = la.lu_factor(L) + T = - la.lu_solve((lu, piv), L_m) + + M_subs = L_0 + np.matmul(L_m, S) + np.matmul(L_p, T) + + return steadystate(Qobj(M_subs, type="super", dims=L_t.dims)) + + def build_preconditioner(A, c_op_list=[], **kwargs): """Constructs a iLU preconditioner necessary for solving for the steady state density matrix using the iterative linear solvers @@ -991,7 +1085,7 @@ def build_preconditioner(A, c_op_list=[], **kwargs): if key in ss_args.keys(): ss_args[key] = kwargs[key] else: - raise Exception("Invalid keyword argument '" + key + + raise TypeError("Invalid keyword argument '" + key + "' passed to steadystate.") # Set column perm to NATURAL if using RCM and not specified by user @@ -1012,7 +1106,7 @@ def build_preconditioner(A, c_op_list=[], **kwargs): ss_list = _steadystate_power_liouvillian(L, ss_args) L, perm, perm2, rev_perm, ss_args = ss_list else: - raise Exception("Invalid preconditioning method.") + raise ValueError("Invalid preconditioning method.") M, ss_args = _iterative_precondition(L, n, ss_args) @@ -1100,7 +1194,7 @@ def _pseudo_inverse_sparse(L, rhoss, w=None, **pseudo_args): A = sp_permute(L.data, perm, perm) Q = sp_permute(Q, perm, perm) else: - if ss_args['solver'] == 'scipy': + if pseudo_args['solver'] == 'scipy': A = L.data.tocsc() A.sort_indices() @@ -1112,10 +1206,9 @@ def _pseudo_inverse_sparse(L, rhoss, w=None, **pseudo_args): else: pspec = pseudo_args['permc_spec'] diag_p_thresh = pseudo_args['diag_pivot_thresh'] - pseudo_args = pseudo_args['ILU_MILU'] lu = sp.linalg.splu(A, permc_spec=pspec, diag_pivot_thresh=diag_p_thresh, - options=dict(ILU_MILU=pseudo_args)) + options=dict(ILU_MILU=pseudo_args['ILU_MILU'])) LIQ = lu.solve(Q.toarray()) elif pseudo_args['method'] == 'spilu': @@ -1136,7 +1229,8 @@ def _pseudo_inverse_sparse(L, rhoss, w=None, **pseudo_args): return Qobj(R, dims=L.dims) -def pseudo_inverse(L, rhoss=None, w=None, sparse=True, **kwargs): +def pseudo_inverse(L, rhoss=None, w=None, sparse=True, + method='splu', **kwargs): """ Compute the pseudo inverse for a Liouvillian superoperator, optionally given its steady state density matrix (which will be computed if not @@ -1194,10 +1288,9 @@ def pseudo_inverse(L, rhoss=None, w=None, sparse=True, **kwargs): if key in pseudo_args.keys(): pseudo_args[key] = kwargs[key] else: - raise Exception( + raise TypeError( "Invalid keyword argument '"+key+"' passed to pseudo_inverse.") - if 'method' not in kwargs.keys(): - pseudo_args['method'] = 'splu' + pseudo_args['method'] = method # Set column perm to NATURAL if using RCM and not specified by user if pseudo_args['use_rcm'] and ('permc_spec' not in kwargs.keys()): diff --git a/qutip/tests/test_control_pulseoptim.py b/qutip/tests/test_control_pulseoptim.py index 69cb4a789f..d71a513429 100644 --- a/qutip/tests/test_control_pulseoptim.py +++ b/qutip/tests/test_control_pulseoptim.py @@ -13,6 +13,7 @@ import tempfile import numpy as np import scipy.optimize +import scipy.sparse as sp import qutip from qutip.control import pulseoptim as cpo @@ -116,6 +117,7 @@ @pytest.fixture(params=[ pytest.param(None, id="default propagation"), pytest.param({'oper_dtype': qutip.Qobj}, id="Qobj propagation"), + pytest.param({'oper_dtype': np.ndarray}, id="ndarray propagation"), ]) def propagation(request): return {'dyn_params': request.param} @@ -165,6 +167,38 @@ def test_basic_optimization(self, system, propagation): result = _optimize_pulse(system) assert result.fid_err < system.kwargs['fid_err_targ'] + def test_sparse_eigen_optimization(self, system, propagation): + system = _merge_kwargs(system, { + **propagation, + "dyn_params": { + **(propagation.get("dyn_params") or {}), + "sparse_eigen_decomp": True, + } + }) + if system.system.dims == [[2], [2]]: + with pytest.raises(ValueError, match=( + r"^Single qubit pulse optimization dynamics cannot use" + r" sparse eigenvector decomposition")): + _optimize_pulse(system) + else: + result = _optimize_pulse(system) + assert result.fid_err < system.kwargs['fid_err_targ'] + + @pytest.mark.parametrize("oper_dtype", [ + pytest.param(sp.csr_matrix, id="csr_matrix"), + pytest.param(list, id="list"), + ]) + def test_invalid_oper_dtype(self, system, oper_dtype): + system = _merge_kwargs(system, {"dyn_params": { + "oper_dtype": oper_dtype, + }}) + with pytest.raises(ValueError) as err: + _optimize_pulse(system) + assert str(err.value) == ( + f"Unknown oper_dtype {oper_dtype!r}. The oper_dtype may be" + " qutip.Qobj or numpy.ndarray." + ) + def test_object_oriented_approach_and_gradient(self, system, propagation): """ Test the object-oriented version of the optimiser, and ensure that the diff --git a/qutip/tests/test_expect.py b/qutip/tests/test_expect.py index b262f08186..4556811071 100644 --- a/qutip/tests/test_expect.py +++ b/qutip/tests/test_expect.py @@ -132,17 +132,29 @@ def test_equivalent_to_matrix_element(hermitian): ]) def test_compatibility_with_solver(solve): e_ops = [getattr(qutip, 'sigma'+x)() for x in 'xyzmp'] + e_ops += [lambda t, psi: np.sin(t)] h = qutip.sigmax() state = qutip.basis(2, 0) times = np.linspace(0, 10, 101) options = qutip.Options(store_states=True) result = solve(h, state, times, e_ops=e_ops, options=options) direct, states = result.expect, result.states - indirect = qutip.expect(e_ops, states) - assert len(direct) == len(indirect) + indirect = qutip.expect(e_ops[:-1], states) + # check measurement operators based on quantum objects + assert len(direct)-1 == len(indirect) for direct_, indirect_ in zip(direct, indirect): assert len(direct_) == len(indirect_) assert isinstance(direct_, np.ndarray) assert isinstance(indirect_, np.ndarray) assert direct_.dtype == indirect_.dtype np.testing.assert_allclose(direct_, indirect_, atol=1e-12) + # test measurement operators based on lambda functions + direct_ = direct[-1] + # by design, lambda measurements are of complex type + indirect_ = np.sin(times, dtype=complex) + assert len(direct_) == len(indirect_) + assert isinstance(direct_, np.ndarray) + assert isinstance(indirect_, np.ndarray) + assert direct_.dtype == indirect_.dtype + np.testing.assert_allclose(direct_, indirect_, atol=1e-12) + diff --git a/qutip/tests/test_gates.py b/qutip/tests/test_gates.py index 588fa27b36..d42ba373e9 100644 --- a/qutip/tests/test_gates.py +++ b/qutip/tests/test_gates.py @@ -299,6 +299,20 @@ def test_cnot_explicit(self): [0, 0, 0, 1, 0, 0, 0, 0]]) np.testing.assert_allclose(test, expected, atol=1e-15) + + def test_hadamard_explicit(self): + test = gates.hadamard_transform(3).full() + expected = np.array([[ 1, 1, 1, 1, 1, 1, 1, 1], + [ 1, -1, 1, -1, 1, -1, 1, -1], + [ 1, 1, -1, -1, 1, 1, -1, -1], + [ 1, -1, -1, 1, 1, -1, -1, 1], + [ 1, 1, 1, 1, -1, -1, -1, -1], + [ 1, -1, 1, -1, -1, 1, -1, 1], + [ 1, 1, -1, -1, -1, -1, 1, 1], + [ 1, -1, -1, 1, -1, 1, 1, -1]]) + expected = expected/np.sqrt(8) + np.testing.assert_allclose(test, expected) + def test_cyclic_permutation(self): operators = [qutip.sigmax(), qutip.sigmaz()] test = gates.expand_operator(qutip.tensor(*operators), N=3, diff --git a/qutip/tests/test_interpolate.py b/qutip/tests/test_interpolate.py index 81f4e5caa3..1028f1abd2 100644 --- a/qutip/tests/test_interpolate.py +++ b/qutip/tests/test_interpolate.py @@ -78,7 +78,7 @@ def _parametrize_solver_coefficients(metafunc): """ size = 10 times = np.linspace(0, 5, 50) - c_ops = [qutip.qzero(size)] + c_ops = [qutip.qeye(size)] solvers = [ (qutip.sesolve, 'sesolve'), (functools.partial(qutip.mesolve, c_ops=c_ops), 'mesolve'), diff --git a/qutip/tests/test_processor.py b/qutip/tests/test_processor.py index 0661b20384..730a165147 100644 --- a/qutip/tests/test_processor.py +++ b/qutip/tests/test_processor.py @@ -142,8 +142,10 @@ def testPlot(self): processor.add_control(sigmaz()) processor.pulses[0].tlist = tlist processor.pulses[0].coeff = np.array([np.sin(t) for t in tlist]) - processor.plot_pulses() - plt.clf() + fig, _ = processor.plot_pulses() + # testing under Xvfb with pytest-xvfb complains if figure windows are + # left open, so we politely close it: + plt.close(fig) # cubic spline tlist = np.linspace(0., 2*np.pi, 20) @@ -151,8 +153,10 @@ def testPlot(self): processor.add_control(sigmaz()) processor.pulses[0].tlist = tlist processor.pulses[0].coeff = np.array([np.sin(t) for t in tlist]) - processor.plot_pulses() - plt.clf() + fig, _ = processor.plot_pulses() + # testing under Xvfb with pytest-xvfb complains if figure windows are + # left open, so we politely close it: + plt.close(fig) def testSpline(self): """ diff --git a/qutip/tests/test_qobj.py b/qutip/tests/test_qobj.py index de20908a96..11caa048bc 100644 --- a/qutip/tests/test_qobj.py +++ b/qutip/tests/test_qobj.py @@ -313,6 +313,11 @@ def test_QobjPower(): q3 = q ** 3 assert (q3.full() - np.linalg.matrix_power(data, 3) < 1e-12).all() +def test_QobjPowerScalar(): + """Check that scalars obtained from bra*ket can be exponentiated. (#1691) + """ + ket = basis(2, 0) + assert (ket.dag()*ket)**2 == Qobj(1) def test_QobjNeg(): "Qobj negation" diff --git a/qutip/tests/test_simdiag.py b/qutip/tests/test_simdiag.py new file mode 100644 index 0000000000..5ece1a9285 --- /dev/null +++ b/qutip/tests/test_simdiag.py @@ -0,0 +1,87 @@ +import pytest +import numpy as np +import qutip + + +@pytest.mark.parametrize('num_mat', [1, 2, 3, 5]) +def test_simdiag(num_mat): + N = 10 + + U = qutip.rand_unitary(N) + commuting_matrices = [U * qutip.qdiags(np.random.rand(N), 0) * U.dag() + for _ in range(num_mat)] + all_evals, evecs = qutip.simdiag(commuting_matrices) + + for matrix, evals in zip(commuting_matrices, all_evals): + for eval, evec in zip(evals, evecs): + assert matrix * evec == evec * eval + + +@pytest.mark.parametrize('num_mat', [1, 2, 3, 5]) +def test_simdiag_no_evals(num_mat): + N = 10 + + U = qutip.rand_unitary(N) + commuting_matrices = [U * qutip.qdiags(np.random.rand(N), 0) * U.dag() + for _ in range(num_mat)] + evecs = qutip.simdiag(commuting_matrices, evals=False) + + for matrix in commuting_matrices: + for evec in evecs: + Mvec = matrix * evec + eval = Mvec.norm() / evec.norm() + assert matrix * evec == evec * eval + + +def test_simdiag_degen(): + N = 10 + U = qutip.rand_unitary(N) + commuting_matrices = [ + U * qutip.qdiags([0, 0, 0, 1, 1, 1, 2, 2, 3, 4], 0) * U.dag(), + U * qutip.qdiags([0, 0, 0, 1, 2, 2, 2, 2, 2, 2], 0) * U.dag(), + U * qutip.qdiags([0, 0, 2, 1, 1, 2, 2, 3, 3, 4], 0) * U.dag(), + ] + all_evals, evecs = qutip.simdiag(commuting_matrices) + + for matrix, evals in zip(commuting_matrices, all_evals): + for eval, evec in zip(evals, evecs): + np.testing.assert_allclose( + (matrix * evec).full(), + (evec * eval).full() + ) + + +@pytest.mark.repeat(10) +def test_simdiag_degen_large(): + N = 20 + U = qutip.rand_unitary(N) + commuting_matrices = [ + U * qutip.qdiags(np.random.randint(0, 3, N), 0) * U.dag() + for _ in range(5) + ] + all_evals, evecs = qutip.simdiag(commuting_matrices, tol=1e-12) + + for matrix, evals in zip(commuting_matrices, all_evals): + for eval, evec in zip(evals, evecs): + np.testing.assert_allclose( + (matrix * evec).full(), + (evec * eval).full() + ) + + +def test_simdiag_no_input(): + with pytest.raises(ValueError): + qutip.simdiag([]) + + +@pytest.mark.parametrize(['ops', 'error'], [ + pytest.param([qutip.basis(5, 0)], 'square', id="Not square"), + pytest.param([qutip.qeye(5), qutip.qeye(3)], 'shape', id="shape mismatch"), + pytest.param([qutip.destroy(5)], 'Hermitian', id="Non Hermitian"), + pytest.param([qutip.sigmax(), qutip.sigmay()], 'commute', + id="Not commuting"), +]) +def test_simdiag_errors(ops, error): + with pytest.raises(TypeError) as err: + qutip.simdiag(ops) + assert error in str(err.value) diff --git a/qutip/tests/test_steadystate.py b/qutip/tests/test_steadystate.py index 941d65f917..549b76d47f 100644 --- a/qutip/tests/test_steadystate.py +++ b/qutip/tests/test_steadystate.py @@ -1,71 +1,39 @@ import numpy as np -from numpy.testing import assert_, assert_equal, run_module_suite - -from qutip import (sigmaz, destroy, steadystate, expect, coherent_dm, - build_preconditioner) - - -def test_qubit_direct(): - "Steady state: Thermal qubit - direct solver" - # thermal steadystate of a qubit: compare numerics with analytical formula - sz = sigmaz() - sm = destroy(2) - - H = 0.5 * 2 * np.pi * sz - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * sm) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * sm.dag()) - rho_ss = steadystate(H, c_op_list, method='direct') - p_ss[idx] = expect(sm.dag() * sm, rho_ss) - - p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec)) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-5, True) - - -def test_qubit_eigen(): - "Steady state: Thermal qubit - eigen solver" - # thermal steadystate of a qubit: compare numerics with analytical formula - sz = sigmaz() - sm = destroy(2) - - H = 0.5 * 2 * np.pi * sz - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * sm) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * sm.dag()) - rho_ss = steadystate(H, c_op_list, method='eigen') - p_ss[idx] = expect(sm.dag() * sm, rho_ss) - - p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec)) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-5, True) - - -def test_qubit_power(): - "Steady state: Thermal qubit - power solver" +import pytest +import qutip +import warnings + + +@pytest.mark.parametrize(['method', 'kwargs'], [ + pytest.param('direct', {}, id="direct"), + pytest.param('direct', {'solver':'mkl'}, id="direct_mkl", + marks=pytest.mark.skipif(not qutip.settings.has_mkl, + reason='MKL extensions not found.')), + pytest.param('direct', {'return_info':True}, id="direct_info"), + pytest.param('direct', {'sparse':False}, id="direct_dense"), + pytest.param('direct', {'use_rcm':True}, id="direct_rcm"), + pytest.param('direct', {'use_wbm':True}, id="direct_wbm"), + pytest.param('eigen', {}, id="eigen"), + pytest.param('eigen', {'use_rcm':True}, id="eigen_rcm"), + pytest.param('svd', {}, id="svd"), + pytest.param('power', {'mtol':1e-5}, id="power"), + pytest.param('power', {'mtol':1e-5, 'solver':'mkl'}, id="power_mkl", + marks=pytest.mark.skipif(not qutip.settings.has_mkl, + reason='MKL extensions not found.')), + pytest.param('power-gmres', {'mtol':1e-1}, id="power-gmres"), + pytest.param('power-gmres', {'mtol':1e-1, 'use_rcm':True, 'use_wbm':True}, + id="power-gmres_perm"), + pytest.param('power-bicgstab', {'use_precond':1}, id="power-bicgstab"), + pytest.param('iterative-gmres', {}, id="iterative-gmres"), + pytest.param('iterative-gmres', {'use_rcm':True, 'use_wbm':True}, + id="iterative-gmres_perm"), + pytest.param('iterative-bicgstab', {'return_info':True}, + id="iterative-bicgstab"), +]) +def test_qubit(method, kwargs): # thermal steadystate of a qubit: compare numerics with analytical formula - sz = sigmaz() - sm = destroy(2) + sz = qutip.sigmaz() + sm = qutip.destroy(2) H = 0.5 * 2 * np.pi * sz gamma1 = 0.05 @@ -73,247 +41,43 @@ def test_qubit_power(): wth_vec = np.linspace(0.1, 3, 20) p_ss = np.zeros(np.shape(wth_vec)) - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * sm) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * sm.dag()) - rho_ss = steadystate(H, c_op_list, method='power', mtol=1e-5) - p_ss[idx] = expect(sm.dag() * sm, rho_ss) + with warnings.catch_warnings(): + if 'use_wbm' in kwargs: + # The deprecation has been fixed in dev.major + warnings.simplefilter("ignore", category=DeprecationWarning) + + for idx, wth in enumerate(wth_vec): + n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature + c_op_list = [] + rate = gamma1 * (1 + n_th) + c_op_list.append(np.sqrt(rate) * sm) + rate = gamma1 * n_th + c_op_list.append(np.sqrt(rate) * sm.dag()) + rho_ss = qutip.steadystate(H, c_op_list, method=method, **kwargs) + if 'return_info' in kwargs: + rho_ss, info = rho_ss + assert isinstance(info, dict) + p_ss[idx] = qutip.expect(sm.dag() * sm, rho_ss) p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec)) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-5, True) - - -def test_qubit_power_gmres(): - "Steady state: Thermal qubit - power-gmres solver" - # thermal steadystate of a qubit: compare numerics with analytical formula - sz = sigmaz() - sm = destroy(2) - - H = 0.5 * 2 * np.pi * sz - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * sm) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * sm.dag()) - rho_ss = steadystate(H, c_op_list, method='power-gmres', mtol=1e-1) - p_ss[idx] = expect(sm.dag() * sm, rho_ss) - - p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec)) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-5, True) - - -def test_qubit_power_bicgstab(): - "Steady state: Thermal qubit - power-bicgstab solver" - # thermal steadystate of a qubit: compare numerics with analytical formula - sz = sigmaz() - sm = destroy(2) - - H = 0.5 * 2 * np.pi * sz - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * sm) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * sm.dag()) - rho_ss = steadystate(H, c_op_list, method='power-bicgstab', use_precond=1) - p_ss[idx] = expect(sm.dag() * sm, rho_ss) - - p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec)) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-5, True) - - - -def test_qubit_gmres(): - "Steady state: Thermal qubit - iterative-gmres solver" - # thermal steadystate of a qubit: compare numerics with analytical formula - sz = sigmaz() - sm = destroy(2) - - H = 0.5 * 2 * np.pi * sz - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * sm) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * sm.dag()) - rho_ss = steadystate(H, c_op_list, method='iterative-gmres') - p_ss[idx] = expect(sm.dag() * sm, rho_ss) - - p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec)) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-5, True) - - -def test_qubit_bicgstab(): - "Steady state: Thermal qubit - iterative-bicgstab solver" - # thermal steadystate of a qubit: compare numerics with analytical formula - sz = sigmaz() - sm = destroy(2) - - H = 0.5 * 2 * np.pi * sz - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * sm) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * sm.dag()) - rho_ss = steadystate(H, c_op_list, method='iterative-bicgstab') - p_ss[idx] = expect(sm.dag() * sm, rho_ss) - - p_ss_analytic = np.exp(-1.0 / wth_vec) / (1 + np.exp(-1.0 / wth_vec)) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-5, True) - - -def test_ho_direct(): - "Steady state: Thermal HO - direct solver" - # thermal steadystate of an oscillator: compare numerics with analytical - # formula - a = destroy(40) - H = 0.5 * 2 * np.pi * a.dag() * a - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * a) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * a.dag()) - rho_ss = steadystate(H, c_op_list, method='direct') - p_ss[idx] = np.real(expect(a.dag() * a, rho_ss)) - - p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-3, True) - - -def test_ho_eigen(): - "Steady state: Thermal HO - eigen solver" - # thermal steadystate of an oscillator: compare numerics with analytical - # formula - a = destroy(40) - H = 0.5 * 2 * np.pi * a.dag() * a - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * a) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * a.dag()) - rho_ss = steadystate(H, c_op_list, method='eigen') - p_ss[idx] = np.real(expect(a.dag() * a, rho_ss)) - - p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-3, True) - - -def test_ho_power(): - "Steady state: Thermal HO - power solver" - # thermal steadystate of an oscillator: compare numerics with analytical - # formula - a = destroy(40) - H = 0.5 * 2 * np.pi * a.dag() * a - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * a) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * a.dag()) - rho_ss = steadystate(H, c_op_list, method='power', mtol=1e-5) - p_ss[idx] = np.real(expect(a.dag() * a, rho_ss)) - - p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-3, True) - -def test_ho_power_gmres(): - "Steady state: Thermal HO - power-gmres solver" - # thermal steadystate of an oscillator: compare numerics with analytical - # formula - a = destroy(40) - H = 0.5 * 2 * np.pi * a.dag() * a - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * a) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * a.dag()) - rho_ss = steadystate(H, c_op_list, method='power-gmres', mtol=1e-1, - use_precond=1) - p_ss[idx] = np.real(expect(a.dag() * a, rho_ss)) - - p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-3, True) - - -def test_ho_power_bicgstab(): - "Steady state: Thermal HO - power-bicgstab solver" + np.testing.assert_allclose(p_ss_analytic, p_ss, atol=1e-5) + + +@pytest.mark.parametrize(['method', 'kwargs'], [ + pytest.param('direct', {}, id="direct"), + pytest.param('direct', {'sparse':False}, id="direct_dense"), + pytest.param('eigen', {}, id="eigen"), + pytest.param('power', {'mtol':1e-5}, id="power"), + pytest.param('power-gmres', {'mtol':1e-1, 'use_precond':1}, id="power-gmres"), + pytest.param('power-bicgstab', {'use_precond':1}, id="power-bicgstab"), + pytest.param('iterative-lgmres', {'use_precond':1}, id="iterative-lgmres"), + pytest.param('iterative-gmres', {}, id="iterative-gmres"), + pytest.param('iterative-bicgstab', {}, id="iterative-bicgstab"), +]) +def test_ho(method, kwargs): # thermal steadystate of an oscillator: compare numerics with analytical # formula - a = destroy(40) + a = qutip.destroy(35) H = 0.5 * 2 * np.pi * a.dag() * a gamma1 = 0.05 @@ -328,188 +92,144 @@ def test_ho_power_bicgstab(): c_op_list.append(np.sqrt(rate) * a) rate = gamma1 * n_th c_op_list.append(np.sqrt(rate) * a.dag()) - rho_ss = steadystate(H, c_op_list, method='power-bicgstab',use_precond=1) - p_ss[idx] = np.real(expect(a.dag() * a, rho_ss)) + rho_ss = qutip.steadystate(H, c_op_list, method=method, **kwargs) + p_ss[idx] = np.real(qutip.expect(a.dag() * a, rho_ss)) p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-3, True) - - -def test_ho_gmres(): - "Steady state: Thermal HO - iterative-gmres solver" - # thermal steadystate of an oscillator: compare numerics with analytical - # formula - a = destroy(40) - H = 0.5 * 2 * np.pi * a.dag() * a - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * a) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * a.dag()) - rho_ss = steadystate(H, c_op_list, method='iterative-gmres') - p_ss[idx] = np.real(expect(a.dag() * a, rho_ss)) - - p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-3, True) - - -def test_ho_bicgstab(): - "Steady state: Thermal HO - iterative-bicgstab solver" - # thermal steadystate of an oscillator: compare numerics with analytical - # formula - a = destroy(40) - H = 0.5 * 2 * np.pi * a.dag() * a - gamma1 = 0.05 - - wth_vec = np.linspace(0.1, 3, 20) - p_ss = np.zeros(np.shape(wth_vec)) - - for idx, wth in enumerate(wth_vec): - - n_th = 1.0 / (np.exp(1.0 / wth) - 1) # bath temperature - c_op_list = [] - rate = gamma1 * (1 + n_th) - c_op_list.append(np.sqrt(rate) * a) - rate = gamma1 * n_th - c_op_list.append(np.sqrt(rate) * a.dag()) - rho_ss = steadystate(H, c_op_list, method='iterative-bicgstab') - p_ss[idx] = np.real(expect(a.dag() * a, rho_ss)) - - p_ss_analytic = 1.0 / (np.exp(1.0 / wth_vec) - 1) - delta = sum(abs(p_ss_analytic - p_ss)) - assert_equal(delta < 1e-3, True) - - - -def test_driven_cavity_direct(): - "Steady state: Driven cavity - direct solver" - - N = 30 - Omega = 0.01 * 2 * np.pi - Gamma = 0.05 - - a = destroy(N) - H = Omega * (a.dag() + a) - c_ops = [np.sqrt(Gamma) * a] - - rho_ss = steadystate(H, c_ops, method='direct') - rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2)) - - assert_((rho_ss - rho_ss_analytic).norm() < 1e-4) - - -def test_driven_cavity_eigen(): - "Steady state: Driven cavity - eigen solver" - + np.testing.assert_allclose(p_ss_analytic, p_ss, atol=1e-3) + + +@pytest.mark.parametrize(['method', 'kwargs'], [ + pytest.param('direct', {}, id="direct"), + pytest.param('direct', {'sparse':False}, id="direct_dense"), + pytest.param('eigen', {}, id="eigen"), + pytest.param('svd', {}, id="svd"), + pytest.param('power', {'mtol':1e-5}, id="power"), + pytest.param('power-gmres', {'mtol':1e-1, 'use_precond':1, 'M':'iterative'}, + id="power-gmres"), + pytest.param('power-bicgstab', {'use_precond':1, 'M':'power'}, + id="power-bicgstab"), + pytest.param('iterative-gmres', {}, id="iterative-gmres"), + pytest.param('iterative-bicgstab', {}, id="iterative-bicgstab"), +]) +def test_driven_cavity(method, kwargs): N = 30 Omega = 0.01 * 2 * np.pi Gamma = 0.05 - a = destroy(N) + a = qutip.destroy(N) H = Omega * (a.dag() + a) c_ops = [np.sqrt(Gamma) * a] - - rho_ss = steadystate(H, c_ops, method='eigen') - rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2)) - - assert_((rho_ss - rho_ss_analytic).norm() < 1e-4) - - -def test_driven_cavity_power(): - "Steady state: Driven cavity - power solver" - - N = 30 - Omega = 0.01 * 2 * np.pi - Gamma = 0.05 - - a = destroy(N) - H = Omega * (a.dag() + a) - c_ops = [np.sqrt(Gamma) * a] - - rho_ss = steadystate(H, c_ops, method='power', mtol=1e-5,) - rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2)) - - assert_((rho_ss - rho_ss_analytic).norm() < 1e-4) - - -def test_driven_cavity_power_gmres(): - "Steady state: Driven cavity - power-gmres solver" - - N = 30 - Omega = 0.01 * 2 * np.pi - Gamma = 0.05 - - a = destroy(N) - H = Omega * (a.dag() + a) - c_ops = [np.sqrt(Gamma) * a] - M = build_preconditioner(H, c_ops, method='power') - rho_ss = steadystate(H, c_ops, method='power-gmres', M=M, mtol=1e-1, - use_precond=1) - rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2)) - assert_((rho_ss - rho_ss_analytic).norm() < 1e-4) - - - -def test_driven_cavity_power_bicgstab(): - "Steady state: Driven cavity - power-bicgstab solver" - - N = 30 - Omega = 0.01 * 2 * np.pi - Gamma = 0.05 - - a = destroy(N) - H = Omega * (a.dag() + a) - c_ops = [np.sqrt(Gamma) * a] - M = build_preconditioner(H, c_ops, method='power') - rho_ss = steadystate(H, c_ops, method='power-bicgstab', M=M, use_precond=1) - rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2)) - assert_((rho_ss - rho_ss_analytic).norm() < 1e-4) - - -def test_driven_cavity_gmres(): - "Steady state: Driven cavity - iterative-gmres solver" - - N = 30 - Omega = 0.01 * 2 * np.pi - Gamma = 0.05 - - a = destroy(N) - H = Omega * (a.dag() + a) - c_ops = [np.sqrt(Gamma) * a] - - rho_ss = steadystate(H, c_ops, method='iterative-gmres') - rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2)) - - assert_((rho_ss - rho_ss_analytic).norm() < 1e-4) - - -def test_driven_cavity_bicgstab(): - "Steady state: Driven cavity - iterative-bicgstab solver" - - N = 30 - Omega = 0.01 * 2 * np.pi - Gamma = 0.05 - - a = destroy(N) - H = Omega * (a.dag() + a) - c_ops = [np.sqrt(Gamma) * a] - - rho_ss = steadystate(H, c_ops, method='iterative-bicgstab') - rho_ss_analytic = coherent_dm(N, -1.0j * (Omega)/(Gamma/2)) - - assert_((rho_ss - rho_ss_analytic).norm() < 1e-4) - - - -if __name__ == "__main__": - run_module_suite() + if 'use_precond' in kwargs: + kwargs['M'] = qutip.build_preconditioner(H, c_ops, method=kwargs['M']) + rho_ss = qutip.steadystate(H, c_ops, method=method, **kwargs) + rho_ss_analytic = qutip.coherent_dm(N, -1.0j * (Omega)/(Gamma/2)) + + np.testing.assert_allclose(rho_ss, rho_ss_analytic, atol=1e-4) + + +@pytest.mark.parametrize(['method', 'kwargs'], [ + pytest.param('splu', {'sparse':False}, id="dense_direct"), + pytest.param('numpy', {'sparse':False}, id="dense_numpy"), + pytest.param('scipy', {'sparse':False}, id="dense_scipy"), + pytest.param('splu', {}, id="splu"), + pytest.param('spilu', {}, id="spilu"), +]) +def test_pseudo_inverse(method, kwargs): + N = 4 + a = qutip.destroy(N) + H = (a.dag() + a) + L = qutip.liouvillian(H, [a]) + rho = qutip.steadystate(L) + Lpinv = qutip.pseudo_inverse(L, rho, method=method, **kwargs) + np.testing.assert_allclose((L * Lpinv * L).full(), L.full()) + np.testing.assert_allclose((Lpinv * L * Lpinv).full(), Lpinv.full()) + + +@pytest.mark.parametrize('sparse', [True, False]) +def test_steadystate_floquet(sparse): + """ + Test the steadystate solution for a periodically + driven system. + """ + N_c = 20 + + a = qutip.destroy(N_c) + a_d = a.dag() + X_c = a + a_d + + w_c = 1 + + A_l = 0.001 + w_l = w_c + gam = 0.01 + + H = w_c * a_d * a + + H_t = [H, [X_c, lambda t, args: args["A_l"] * np.cos(args["w_l"] * t)]] + + psi0 = qutip.fock(N_c, 0) + + args = {"A_l": A_l, "w_l": w_l} + + c_ops = [] + c_ops.append(np.sqrt(gam) * a) + + t_l = np.linspace(0, 20 / gam, 2000) + + expect_me = qutip.mesolve(H_t, psi0, t_l, + c_ops, [a_d * a], args=args).expect[0] + + rho_ss = qutip.steadystate_floquet(H, c_ops, + A_l * X_c, w_l, n_it=3, sparse=sparse) + expect_ss = qutip.expect(a_d * a, rho_ss) + + np.testing.assert_allclose(expect_me[-20:], expect_ss, atol=1e-3) + + +def test_bad_options_steadystate(): + N = 4 + a = qutip.destroy(N) + H = (a.dag() + a) + c_ops = [a] + with pytest.raises(ValueError): + rho_ss = qutip.steadystate(H, c_ops, method='not a method') + with pytest.raises(TypeError): + rho_ss = qutip.steadystate(H, c_ops, method='direct', bad_opt=True) + with pytest.raises(ValueError): + rho_ss = qutip.steadystate(H, c_ops, method='direct', solver='Error') + + +def test_bad_options_pseudo_inverse(): + N = 4 + a = qutip.destroy(N) + H = (a.dag() + a) + L = qutip.liouvillian(H, [a]) + with pytest.raises(TypeError): + qutip.pseudo_inverse(L, method='splu', bad_opt=True) + with pytest.raises(ValueError): + qutip.pseudo_inverse(L, method='not a method', sparse=False) + with pytest.raises(ValueError): + qutip.pseudo_inverse(L, method='not a method') + + +def test_bad_options_build_preconditioner(): + N = 4 + a = qutip.destroy(N) + H = (a.dag() + a) + c_ops = [a] + with pytest.raises(TypeError): + qutip.build_preconditioner(H, c_ops, method='power', bad_opt=True) + with pytest.raises(ValueError): + qutip.build_preconditioner(H, c_ops, method='not a method') + + +def test_bad_system(): + N = 4 + a = qutip.destroy(N) + H = (a.dag() + a) + c_ops = [a] + with pytest.raises(TypeError) as err: + rho_ss = qutip.steadystate(H, [], method='direct') + with pytest.raises(TypeError) as err: + rho_ss = qutip.steadystate(qutip.basis(N, N-1), [], method='direct') diff --git a/qutip/tests/test_utilities.py b/qutip/tests/test_utilities.py index 9961f2018c..e40c3b6800 100644 --- a/qutip/tests/test_utilities.py +++ b/qutip/tests/test_utilities.py @@ -1,104 +1,107 @@ import numpy as np -from numpy.testing import assert_, run_module_suite +from qutip import convert_unit, clebsch, n_thermal +import qutip.utilities as utils +from functools import partial +import pytest -from qutip import convert_unit, clebsch +@pytest.mark.parametrize(['w', 'w_th', 'expected'], [ + pytest.param(np.log(2), 1, 1, id='log(2)'), + pytest.param(np.log(2)*5, 5, 1, id='5*log(2)'), + pytest.param(0, 1, 0, id="0_energy"), + pytest.param(1, -1, 0, id="neg_temp"), + pytest.param(np.array([np.log(2), np.log(3), np.log(4)]), 1, + np.array([1, 1/2, 1/3]), id="array"), +]) +def test_n_thermal(w, w_th, expected): + np.testing.assert_allclose(n_thermal(w, w_th), expected) -def test_unit_conversions(): - "utilities: energy unit conversions" +def _get_converter(orig, target): + """get funtion 'convert_{}_to_{}' when available for coverage """ + try: + func = getattr(utils, f'convert_{orig}_to_{target}') + except AttributeError: + func = partial(convert_unit, orig=orig, to=target) + return func + + +@pytest.mark.parametrize('orig', ["J", "eV", "meV", "GHz", "mK"]) +@pytest.mark.parametrize('target', ["J", "eV", "meV", "GHz", "mK"]) +def test_unit_conversions(orig, target): + T = np.random.rand() * 100.0 + T_converted = convert_unit(T, orig=orig, to=target) + T_back = convert_unit(T_converted, orig=target, to=orig) + + assert T == pytest.approx(T_back) + + T_converted = _get_converter(orig=orig, target=target)(T) + T_back = _get_converter(orig=target, target=orig)(T_converted) + + assert T == pytest.approx(T_back) + + +@pytest.mark.parametrize('orig', ["J", "eV", "meV", "GHz", "mK"]) +@pytest.mark.parametrize('middle', ["J", "eV", "meV", "GHz", "mK"]) +@pytest.mark.parametrize('target', ["J", "eV", "meV", "GHz", "mK"]) +def test_unit_conversions_loop(orig, middle, target): T = np.random.rand() * 100.0 + T_middle = convert_unit(T, orig=orig, to=middle) + T_converted = convert_unit(T_middle, orig=middle, to=target) + T_back = convert_unit(T_converted, orig=target, to=orig) + + assert T == pytest.approx(T_back) + - diff = convert_unit(convert_unit(T, orig="mK", to="GHz"), - orig="GHz", to="mK") - T - assert_(abs(diff) < 1e-6) - diff = convert_unit(convert_unit(T, orig="mK", to="meV"), - orig="meV", to="mK") - T - assert_(abs(diff) < 1e-6) - - diff = convert_unit(convert_unit(convert_unit(T, orig="mK", to="GHz"), - orig="GHz", to="meV"), - orig="meV", to="mK") - T - assert_(abs(diff) < 1e-6) - - w = np.random.rand() * 100.0 - - diff = convert_unit(convert_unit(w, orig="GHz", to="meV"), - orig="meV", to="GHz") - w - assert_(abs(diff) < 1e-6) - - diff = convert_unit(convert_unit(w, orig="GHz", to="mK"), - orig="mK", to="GHz") - w - assert_(abs(diff) < 1e-6) - - diff = convert_unit(convert_unit(convert_unit(w, orig="GHz", to="mK"), - orig="mK", to="meV"), - orig="meV", to="GHz") - w - assert_(abs(diff) < 1e-6) - -def test_unit_clebsch(): - "utilities: Clebsch–Gordan coefficients " - N = 15 - for _ in range(100): - "sum_m1 sum_m2 C(j1,j2,j3,m1,m2,m3)*C(j1,j2,j3',m1,m2,m3') =" - "delta j3,j3' delta m3,m3'" - j1 = np.random.randint(0, N+1) - j2 = np.random.randint(0, N+1) - j3 = np.random.randint(abs(j1-j2), j1+j2+1) - j3p = np.random.randint(abs(j1-j2), j1+j2+1) - m3 = np.random.randint(-j3, j3+1) - m3p = np.random.randint(-j3p, j3p+1) - if np.random.rand() < 0.25: - j1 += 0.5 - j3 += 0.5 - j3p += 0.5 - m3 += np.random.choice([-0.5, 0.5]) - m3p += np.random.choice([-0.5, 0.5]) - if np.random.rand() < 0.25: - j2 += 0.5 - j3 += 0.5 - j3p += 0.5 - m3 += np.random.choice([-0.5, 0.5]) - m3p += np.random.choice([-0.5, 0.5]) - sum_match = -1 - sum_differ = -int(j3 == j3p and m3 == m3p) - for m1 in np.arange(-j1,j1+1): - for m2 in np.arange(-j2,j2+1): +def test_unit_conversions_bad_unit(): + with pytest.raises(TypeError): + convert_unit(10, orig="bad", to="J") + with pytest.raises(TypeError): + convert_unit(10, orig="J", to="bad") + + + +@pytest.mark.parametrize('j1', [0.5, 1.0, 1.5, 2.0, 5, 7.5, 10, 12.5]) +@pytest.mark.parametrize('j2', [0.5, 1.0, 1.5, 2.0, 5, 7.5, 10, 12.5]) +def test_unit_clebsch_delta_j(j1, j2): + """sum_m1 sum_m2 C(j1,j2,j3,m1,m2,m3) * C(j1,j2,j3',m1,m2,m3') = + delta j3,j3' delta m3,m3'""" + for _ in range(10): + j3 = np.random.choice(np.arange(abs(j1-j2), j1+j2+1)) + j3p = np.random.choice(np.arange(abs(j1-j2), j1+j2+1)) + m3 = np.random.choice(np.arange(-j3, j3+1)) + m3p = np.random.choice(np.arange(-j3p, j3p+1)) + + sum_match = 0 + sum_differ = 0 + for m1 in np.arange(-j1, j1+1): + for m2 in np.arange(-j2, j2+1): c1 = clebsch(j1, j2, j3, m1, m2, m3) c2 = clebsch(j1, j2, j3p, m1, m2, m3p) sum_match += c1**2 sum_differ += c1*c2 - assert_(abs(sum_match) < 1e-6) - assert_(abs(sum_differ) < 1e-6) - - for _ in range(100): - "sum_j3 sum_m3 C(j1,j2,j3,m1,m2,m3)*C(j1,j2,j3,m1',m2',m3) =" - "delta m1,m1' delta m2,m2'" - j1 = np.random.randint(0,N+1) - j2 = np.random.randint(0,N+1) - m1 = np.random.randint(-j1,j1+1) - m1p = np.random.randint(-j1,j1+1) - m2 = np.random.randint(-j2,j2+1) - m2p = np.random.randint(-j2,j2+1) - if np.random.rand() < 0.25: - j1 += 0.5 - m1 += np.random.choice([-0.5, 0.5]) - m1p += np.random.choice([-0.5, 0.5]) - if np.random.rand() < 0.25: - j2 += 0.5 - m2 += np.random.choice([-0.5, 0.5]) - m2p += np.random.choice([-0.5, 0.5]) - sum_match = -1 - sum_differ = -int(m1 == m1p and m2 == m2p) - for j3 in np.arange(abs(j1-j2),j1+j2+1): - for m3 in np.arange(-j3,j3+1): + assert sum_match == pytest.approx(1) + assert sum_differ == pytest.approx(int(j3 == j3p and m3 == m3p)) + + +@pytest.mark.parametrize('j1', [0.5, 1.0, 1.5, 2.0, 5, 7.5, 10, 12.5]) +@pytest.mark.parametrize('j2', [0.5, 1.0, 1.5, 2.0, 5, 7.5, 10, 12.5]) +def test_unit_clebsch_delta_m(j1, j2): + """sum_j3 sum_m3 C(j1,j2,j3,m1,m2,m3)*C(j1,j2,j3,m1',m2',m3) = + delta m1,m1' delta m2,m2'""" + for _ in range(10): + m1 = np.random.choice(np.arange(-j1, j1+1)) + m1p = np.random.choice(np.arange(-j1, j1+1)) + m2 = np.random.choice(np.arange(-j2, j2+1)) + m2p = np.random.choice(np.arange(-j2, j2+1)) + + sum_match = 0 + sum_differ = 0 + for j3 in np.arange(abs(j1-j2), j1+j2+1): + for m3 in np.arange(-j3, j3+1): c1 = clebsch(j1, j2, j3, m1, m2, m3) c2 = clebsch(j1, j2, j3, m1p, m2p, m3) sum_match += c1**2 sum_differ += c1*c2 - assert_(abs(sum_match) < 1e-6) - assert_(abs(sum_differ) < 1e-6) - - -if __name__ == "__main__": - run_module_suite() + assert sum_match == pytest.approx(1) + assert sum_differ == pytest.approx(int(m1 == m1p and m2 == m2p)) diff --git a/qutip/tests/test_wigner.py b/qutip/tests/test_wigner.py index 7c1c63bc5c..26c2030358 100644 --- a/qutip/tests/test_wigner.py +++ b/qutip/tests/test_wigner.py @@ -1,8 +1,10 @@ import pytest import numpy as np +import itertools from scipy.special import laguerre from numpy.random import rand -from numpy.testing import assert_, run_module_suite, assert_equal +from numpy.testing import assert_, run_module_suite, assert_equal, \ + assert_almost_equal, assert_allclose import qutip from qutip.states import coherent, fock, ket, bell_state @@ -593,6 +595,102 @@ def test_wigner_clenshaw_sp_iter_dm(): Wdiff = abs(W - Wclen) assert_equal(np.sum(abs(Wdiff)) < 1e-7, True) +@pytest.mark.parametrize(['spin'], [ + pytest.param(1/2, id="spin-one-half"), + pytest.param(3, id="spin-three"), + pytest.param(13/2, id="spin-thirteen-half"), + pytest.param(7, id="spin-seven") +]) +@pytest.mark.parametrize("pure", [ + pytest.param(True, id="pure"), + pytest.param(False, id="mixed") +]) +def test_spin_q_function(spin, pure): + d = int(2*spin + 1) + rho = rand_dm(d, pure=pure) + + # Points at which to evaluate the spin Q function + theta = np.linspace(0, np.pi, 16, endpoint=True) + phi = np.linspace(-np.pi, np.pi, 32, endpoint=True) + Q, _, _ = qutip.spin_q_function(rho, theta, phi) + + for k, (phi_prime, theta_prime) in enumerate(itertools.product(phi, theta)): + state = qutip.spin_coherent(spin, theta_prime, phi_prime) + direct_Q = (state.dag() * rho * state).norm() + assert_almost_equal(Q.flat[k], direct_Q, decimal=9) + +@pytest.mark.parametrize(['spin'], [ + pytest.param(1/2, id="spin-one-half"), + pytest.param(3, id="spin-three"), + pytest.param(13/2, id="spin-thirteen-half"), + pytest.param(7, id="spin-seven") +]) +@pytest.mark.parametrize("pure", [ + pytest.param(True, id="pure"), + pytest.param(False, id="mixed") +]) +def test_spin_q_function_normalized(spin, pure): + d = int(2 * spin + 1) + rho = rand_dm(d, pure=pure) + + # Points at which to evaluate the spin Q function + theta = np.linspace(0, np.pi, 128, endpoint=True) + phi = np.linspace(-np.pi, np.pi, 256, endpoint=True) + Q, THETA, _ = qutip.spin_q_function(rho, theta, phi) + + norm = d / (4 * np.pi) * np.trapz(np.trapz(Q * np.sin(THETA), theta), phi) + assert_almost_equal(norm, 1, decimal=4) + + +@pytest.mark.parametrize(["spin"], [ + pytest.param(1/2, id="spin-one-half"), + pytest.param(1, id="spin-one"), + pytest.param(3/2, id="spin-three-half"), + pytest.param(2, id="spin-two") +]) +@pytest.mark.parametrize("pure", [ + pytest.param(True, id="pure"), + pytest.param(False, id="mixed") +]) +def test_spin_wigner_normalized(spin, pure): + d = int(2*spin + 1) + rho = rand_dm(d, pure=pure) + + # Points at which to evaluate the spin Wigner function + theta = np.linspace(0, np.pi, 256, endpoint=True) + phi = np.linspace(-np.pi, np.pi, 512, endpoint=True) + W, THETA, PHI = qutip.spin_wigner(rho, theta, phi) + + norm = np.trapz(np.trapz(W * np.sin(THETA) * np.sqrt(d / (4*np.pi)), theta), phi) + assert_almost_equal(norm, 1, decimal=4) + +@pytest.mark.parametrize(['spin'], [ + pytest.param(1 / 2, id="spin-one-half"), + pytest.param(1, id="spin-one"), + pytest.param(3 / 2, id="spin-three-half"), + pytest.param(2, id="spin-two") +]) +@pytest.mark.parametrize("pure", [ + pytest.param(True, id="pure"), + pytest.param(False, id="mixed") +]) +def test_spin_wigner_overlap(spin, pure, n=5): + d = int(2*spin + 1) + rho = rand_dm(d, pure=pure) + + # Points at which to evaluate the spin Wigner function + theta = np.linspace(0, np.pi, 256, endpoint=True) + phi = np.linspace(-np.pi, np.pi, 512, endpoint=True) + W, THETA, _ = qutip.spin_wigner(rho, theta, phi) + + for k in range(n): + test_state = rand_dm(d) + state_overlap = (test_state*rho).tr().real + + W_state, _, _ = qutip.spin_wigner(test_state, theta, phi) + W_overlap = np.trapz( + np.trapz(W_state * W * np.sin(THETA), theta), phi).real + assert_almost_equal(W_overlap, state_overlap, decimal=4) if __name__ == "__main__": run_module_suite() diff --git a/qutip/utilities.py b/qutip/utilities.py index 5740961fd9..ba7f5c8d69 100644 --- a/qutip/utilities.py +++ b/qutip/utilities.py @@ -3,8 +3,7 @@ qutip modules. """ -__all__ = ['n_thermal', 'linspace_with', 'clebsch', 'convert_unit', - 'view_methods'] +__all__ = ['n_thermal', 'clebsch', 'convert_unit'] import numpy as np @@ -45,39 +44,6 @@ def n_thermal(w, w_th): return 0.0 -def linspace_with(start, stop, num=50, elems=[]): - """ - Return an array of numbers sampled over specified interval - with additional elements added. - - Returns `num` spaced array with elements from `elems` inserted - if not already included in set. - - Returned sample array is not evenly spaced if addtional elements - are added. - - Parameters - ---------- - start : int - The starting value of the sequence. - stop : int - The stoping values of the sequence. - num : int, optional - Number of samples to generate. - elems : list/ndarray, optional - Requested elements to include in array - - Returns - ------- - samples : ndadrray - Original equally spaced sample array with additional - elements added. - """ - elems = np.array(elems) - lspace = np.linspace(start, stop, num) - return np.union1d(lspace, elems) - - def _factorial_prod(N, arr): arr[:int(N)] += 1 @@ -358,32 +324,6 @@ def convert_mK_to_GHz(w): return w_GHz -def view_methods(Q): - """ - View the methods and corresponding doc strings - for a Qobj class. - - Parameters - ---------- - Q : Qobj - Input Quantum object. - - """ - meth = dir(Q) - qobj_props = ['data', 'dims', 'isherm', 'shape', 'type'] - pub_meth = [x for x in meth if x.find('_') and x not in qobj_props] - ml = max([len(x) for x in pub_meth]) - nl = len(Q.__class__.__name__ + 'Class Methods:') - print(Q.__class__.__name__ + ' Class Methods:') - print('-' * nl) - for ii in range(len(pub_meth)): - m = getattr(Q, pub_meth[ii]) - meth_str = m.__doc__ - ind = meth_str.find('\n') - pub_len = len(pub_meth[ii] + ': ') - print(pub_meth[ii] + ':' + ' ' * (ml+3-pub_len) + meth_str[:ind]) - - def _version2int(version_string): str_list = version_string.split( "-dev")[0].split("rc")[0].split("a")[0].split("b")[0].split( diff --git a/qutip/visualization.py b/qutip/visualization.py index 616e4318e7..47f40fc101 100644 --- a/qutip/visualization.py +++ b/qutip/visualization.py @@ -382,8 +382,103 @@ def sphereplot(theta, phi, values, fig=None, ax=None, save=False): return fig, ax +def _remove_margins(axis): + """ + removes margins about z = 0 and improves the style + by monkey patching + """ + def _get_coord_info_new(renderer): + mins, maxs, centers, deltas, tc, highs = \ + _get_coord_info_old(renderer) + mins += deltas / 4 + maxs -= deltas / 4 + return mins, maxs, centers, deltas, tc, highs + + _get_coord_info_old = axis._get_coord_info + axis._get_coord_info = _get_coord_info_new + + +def _truncate_colormap(cmap, minval=0.0, maxval=1.0, n=100): + """ + truncates portion of a colormap and returns the new one + """ + if isinstance(cmap, str): + cmap = plt.get_cmap(cmap) + new_cmap = mpl.colors.LinearSegmentedColormap.from_list( + 'trunc({n},{a:.2f},{b:.2f})'.format( + n=cmap.name, a=minval, b=maxval), + cmap(np.linspace(minval, maxval, n))) + return new_cmap + + +def _stick_to_planes(stick, azim, ax, M, spacing): + """adjusts xlim and ylim in way that bars will + Stick to xz and yz planes + """ + if stick is True: + azim = azim % 360 + if 0 <= azim <= 90: + ax.set_ylim(1 - .5,) + ax.set_xlim(1 - .5,) + elif 90 < azim <= 180: + ax.set_ylim(1 - .5,) + ax.set_xlim(0, M.shape[0] + (.5 - spacing)) + elif 180 < azim <= 270: + ax.set_ylim(0, M.shape[1] + (.5 - spacing)) + ax.set_xlim(0, M.shape[0] + (.5 - spacing)) + elif 270 < azim < 360: + ax.set_ylim(0, M.shape[1] + (.5 - spacing)) + ax.set_xlim(1 - .5,) + + +def _update_yaxis(spacing, M, ax, ylabels): + """ + updates the y-axis + """ + ytics = [x + (1 - (spacing / 2)) for x in range(M.shape[1])] + ax.axes.w_yaxis.set_major_locator(plt.FixedLocator(ytics)) + if ylabels: + nylabels = len(ylabels) + if nylabels != len(ytics): + raise ValueError(f"got {nylabels} ylabels but needed {len(ytics)}") + ax.set_yticklabels(ylabels) + else: + ax.set_yticklabels([str(y + 1) for y in range(M.shape[1])]) + ax.set_yticklabels([str(i) for i in range(M.shape[1])]) + ax.tick_params(axis='y', labelsize=14) + ax.set_yticks([y + (1 - (spacing / 2)) for y in range(M.shape[1])]) + + +def _update_xaxis(spacing, M, ax, xlabels): + """ + updates the x-axis + """ + xtics = [x + (1 - (spacing / 2)) for x in range(M.shape[1])] + ax.axes.w_xaxis.set_major_locator(plt.FixedLocator(xtics)) + if xlabels: + nxlabels = len(xlabels) + if nxlabels != len(xtics): + raise ValueError(f"got {nxlabels} xlabels but needed {len(xtics)}") + ax.set_xticklabels(xlabels) + else: + ax.set_xticklabels([str(x + 1) for x in range(M.shape[0])]) + ax.set_xticklabels([str(i) for i in range(M.shape[0])]) + ax.tick_params(axis='x', labelsize=14) + ax.set_xticks([x + (1 - (spacing / 2)) for x in range(M.shape[0])]) + + +def _update_zaxis(ax, z_min, z_max, zticks): + """ + updates the z-axis + """ + ax.axes.w_zaxis.set_major_locator(plt.IndexLocator(1, 0.5)) + if isinstance(zticks, list): + ax.set_zticks(zticks) + ax.set_zlim3d([min(z_min, 0), z_max]) + + def matrix_histogram(M, xlabels=None, ylabels=None, title=None, limits=None, - colorbar=True, fig=None, ax=None): + colorbar=True, fig=None, ax=None, options=None): """ Draw a histogram for the matrix M, with the given x and y labels and title. @@ -407,7 +502,74 @@ def matrix_histogram(M, xlabels=None, ylabels=None, title=None, limits=None, ax : a matplotlib axes instance The axes context in which the plot will be drawn. - Returns + colorbar : bool (default: True) + show colorbar + + options : dict + A dictionary containing extra options for the plot. + The names (keys) and values of the options are + described below: + + 'zticks' : list of numbers + A list of z-axis tick locations. + + 'cmap' : string (default: 'jet') + The name of the color map to use. + + 'cmap_min' : float (default: 0.0) + The lower bound to truncate the color map at. + A value in range 0 - 1. The default, 0, leaves the lower + bound of the map unchanged. + + 'cmap_max' : float (default: 1.0) + The upper bound to truncate the color map at. + A value in range 0 - 1. The default, 1, leaves the upper + bound of the map unchanged. + + 'bars_spacing' : float (default: 0.1) + spacing between bars. + + 'bars_alpha' : float (default: 1.) + transparency of bars, should be in range 0 - 1 + + 'bars_lw' : float (default: 0.5) + linewidth of bars' edges. + + 'bars_edgecolor' : color (default: 'k') + The colors of the bars' edges. + Examples: 'k', (0.1, 0.2, 0.5) or '#0f0f0f80'. + + 'shade' : bool (default: True) + Whether to shade the dark sides of the bars (True) or not (False). + The shading is relative to plot's source of light. + + 'azim' : float + The azimuthal viewing angle. + + 'elev' : float + The elevation viewing angle. + + 'proj_type' : string (default: 'ortho' if ax is not passed) + The type of projection ('ortho' or 'persp') + + 'stick' : bool (default: False) + Changes xlim and ylim in such a way that bars next to + XZ and YZ planes will stick to those planes. + This option has no effect if ``ax`` is passed as a parameter. + + 'cbar_pad' : float (default: 0.04) + The fraction of the original axes between the colorbar + and the new image axes. + (i.e. the padding between the 3D figure and the colorbar). + + 'cbar_to_z' : bool (default: False) + Whether to set the color of maximum and minimum z-values to the + maximum and minimum colors in the colorbar (True) or not (False). + + 'figsize' : tuple of two numbers + The size of the figure. + + Returns : ------- fig, ax : tuple A tuple of the matplotlib figure and axes instances used to produce @@ -420,16 +582,38 @@ def matrix_histogram(M, xlabels=None, ylabels=None, title=None, limits=None, """ + # default options + default_opts = {'figsize': None, 'cmap': 'jet', 'cmap_min': 0., + 'cmap_max': 1., 'zticks': None, 'bars_spacing': 0.2, + 'bars_alpha': 1., 'bars_lw': 0.5, 'bars_edgecolor': 'k', + 'shade': False, 'azim': -35, 'elev': 35, + 'proj_type': 'ortho', 'stick': False, + 'cbar_pad': 0.04, 'cbar_to_z': False} + + # update default_opts from input options + if options is None: + pass + elif isinstance(options, dict): + # check if keys in options dict are valid + if set(options) - set(default_opts): + raise ValueError("invalid key(s) found in options: " + f"{', '.join(set(options) - set(default_opts))}") + else: + # updating default options + default_opts.update(options) + else: + raise ValueError("options must be a dictionary") + if isinstance(M, Qobj): # extract matrix data from Qobj M = M.full() n = np.size(M) xpos, ypos = np.meshgrid(range(M.shape[0]), range(M.shape[1])) - xpos = xpos.T.flatten() - 0.5 - ypos = ypos.T.flatten() - 0.5 + xpos = xpos.T.flatten() + 0.5 + ypos = ypos.T.flatten() + 0.5 zpos = np.zeros(n) - dx = dy = 0.8 * np.ones(n) + dx = dy = (1 - default_opts['bars_spacing']) * np.ones(n) dz = np.real(M.flatten()) if isinstance(limits, list) and len(limits) == 2: @@ -442,48 +626,59 @@ def matrix_histogram(M, xlabels=None, ylabels=None, title=None, limits=None, z_min -= 0.1 z_max += 0.1 - norm = mpl.colors.Normalize(z_min, z_max) - cmap = cm.get_cmap('jet') # Spectral + if default_opts['cbar_to_z']: + norm = mpl.colors.Normalize(min(dz), max(dz)) + else: + norm = mpl.colors.Normalize(z_min, z_max) + cmap = _truncate_colormap(default_opts['cmap'], + default_opts['cmap_min'], + default_opts['cmap_max']) colors = cmap(norm(dz)) if ax is None: - fig = plt.figure() - ax = _axes3D(fig, azim=-35, elev=35) - - ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color=colors) + fig = plt.figure(figsize=default_opts['figsize']) + ax = _axes3D(fig, + azim=default_opts['azim'] % 360, + elev=default_opts['elev'] % 360) + ax.set_proj_type(default_opts['proj_type']) + + ax.bar3d(xpos, ypos, zpos, dx, dy, dz, color=colors, + edgecolors=default_opts['bars_edgecolor'], + linewidths=default_opts['bars_lw'], + alpha=default_opts['bars_alpha'], + shade=default_opts['shade']) + # remove vertical lines on xz and yz plane + ax.yaxis._axinfo["grid"]['linewidth'] = 0 + ax.xaxis._axinfo["grid"]['linewidth'] = 0 if title and fig: ax.set_title(title) # x axis - xtics = -0.5 + np.arange(M.shape[0]) - ax.axes.w_xaxis.set_major_locator(plt.FixedLocator(xtics)) - if xlabels: - nxlabels = len(xlabels) - if nxlabels != len(xtics): - raise ValueError(f"got {nxlabels} xlabels but needed {len(xtics)}") - ax.set_xticklabels(xlabels) - ax.tick_params(axis='x', labelsize=14) + _update_xaxis(default_opts['bars_spacing'], M, ax, xlabels) # y axis - ytics = -0.5 + np.arange(M.shape[1]) - ax.axes.w_yaxis.set_major_locator(plt.FixedLocator(ytics)) - if ylabels: - nylabels = len(ylabels) - if nylabels != len(ytics): - raise ValueError(f"got {nylabels} ylabels but needed {len(ytics)}") - ax.set_yticklabels(ylabels) - ax.tick_params(axis='y', labelsize=14) + _update_yaxis(default_opts['bars_spacing'], M, ax, ylabels) # z axis - ax.axes.w_zaxis.set_major_locator(plt.IndexLocator(1, 0.5)) - ax.set_zlim3d([min(z_min, 0), z_max]) + _update_zaxis(ax, z_min, z_max, default_opts['zticks']) + + # stick to xz and yz plane + _stick_to_planes(default_opts['stick'], + default_opts['azim'], ax, M, + default_opts['bars_spacing']) # color axis if colorbar: - cax, kw = mpl.colorbar.make_axes(ax, shrink=.75, pad=.0) + cax, kw = mpl.colorbar.make_axes(ax, shrink=.75, + pad=default_opts['cbar_pad']) mpl.colorbar.ColorbarBase(cax, cmap=cmap, norm=norm) + # removing margins + _remove_margins(ax.xaxis) + _remove_margins(ax.yaxis) + _remove_margins(ax.zaxis) + return fig, ax diff --git a/qutip/wigner.py b/qutip/wigner.py index 6b85a47708..abe987fd3a 100644 --- a/qutip/wigner.py +++ b/qutip/wigner.py @@ -850,14 +850,24 @@ def qfunc( # PSEUDO DISTRIBUTION FUNCTIONS FOR SPINS # def spin_q_function(rho, theta, phi): - """Husimi Q-function for spins. + r"""The Husimi Q function for spins is defined as ``Q(theta, phi) = + SCS.dag() * rho * SCS`` for the spin coherent state ``SCS = spin_coherent( + j, theta, phi)`` where j is the spin length. + The implementation here is more efficient as it doesn't + generate all of the SCS at theta and phi (see references). + + The spin Q function is normal when integrated over the surface of the + sphere + + .. math:: \frac{4 \pi}{2j + 1}\int_\phi \int_\theta + Q(\theta, \phi) \sin(\theta) d\theta d\phi = 1 Parameters ---------- state : qobj A state vector or density matrix for a spin-j quantum system. theta : array_like - Polar angle at which to calculate the Husimi-Q function. + Polar (colatitude) angle at which to calculate the Husimi-Q function. phi : array_like Azimuthal angle at which to calculate the Husimi-Q function. @@ -867,6 +877,11 @@ def spin_q_function(rho, theta, phi): Values representing the spin Husimi Q function at the values specified by THETA and PHI. + References + ---------- + [1] Lee Loh, Y., & Kim, M. (2015). American J. of Phys., 83(1), 30–35. + https://doi.org/10.1119/1.4898595 + """ if rho.type == 'bra': @@ -882,42 +897,71 @@ def spin_q_function(rho, theta, phi): Q = np.zeros_like(THETA, dtype=complex) - for m1 in arange(-j, j+1): - - Q += binom(2*j, j+m1) * cos(THETA/2) ** (2*(j-m1)) * sin(THETA/2) ** (2*(j+m1)) * \ - rho.data[int(j-m1), int(j-m1)] + for m1 in arange(-j, j + 1): + Q += binom(2 * j, j + m1) * cos(THETA / 2) ** (2 * (j + m1)) * \ + sin(THETA / 2) ** (2 * (j - m1)) * \ + rho.data[int(j - m1), int(j - m1)] - for m2 in arange(m1+1, j+1): + for m2 in arange(m1 + 1, j + 1): + Q += (sqrt(binom(2 * j, j + m1)) * sqrt(binom(2 * j, j + m2)) * + cos(THETA / 2) ** (2 * j + m1 + m2) * + sin(THETA / 2) ** (2 * j - m1 - m2)) * \ + (exp(1j * (m1 - m2) * PHI) * rho.data[int(j - m1), int(j - m2)] + + exp(1j * (m2 - m1) * PHI) * rho.data[int(j - m2), int(j - m1)]) - Q += (sqrt(binom(2*j, j+m1)) * sqrt(binom(2*j, j+m2)) * - cos(THETA/2) ** (2*j-m1-m2) * sin(THETA/2) ** (2*j+m1+m2)) * \ - (exp(1j * (m2-m1) * PHI) * rho.data[int(j-m1), int(j-m2)] + - exp(1j * (m1-m2) * PHI) * rho.data[int(j-m2), int(j-m1)]) + return Q.real, THETA, PHI - return Q.real/pi, THETA, PHI def _rho_kq(rho, j, k, q): + """ + This calculates the trace of the multipole operator T_kq and the density + matrix rho for use in the spin Wigner quasiprobability distribution. + + Parameters + ---------- + rho : qobj + A density matrix for a spin-j quantum system. + j : float + The spin length of the system. + k : int + Spherical harmonic degree + q : int + Spherical harmonic order + + Returns + ------- + v : float + Overlap of state with multipole operator T_kq + """ + v = 0j for m1 in arange(-j, j+1): for m2 in arange(-j, j+1): v += ( - (-1)**(j - m1 - q) - * qutip.clebsch(j, j, k, m1, -m2, q) - * rho.data[m1 + j, m2 + j] + (-1) ** (2 * j - k - m1 - m2) + * np.sqrt((2 * k + 1) / (2 * j + 1)) + * qutip.clebsch(j, k, j, -m1, q, -m2) + * rho.data[int(j - m1), int(j - m2)] ) return v def spin_wigner(rho, theta, phi): - """Wigner function for a spin-j system on the 2-sphere of radius j - (for j = 1/2 this is the Bloch sphere). + r"""Wigner function for a spin-j system. + + The spin W function is normal when integrated over the surface of the + sphere + + .. math:: \sqrt{\frac{4 \pi}{2j + 1}}\int_\phi \int_\theta + W(\theta,\phi) \sin(\theta) d\theta d\phi = 1 + Parameters ---------- state : qobj A state vector or density matrix for a spin-j quantum system. theta : array_like - Polar angle at which to calculate the W function. + Polar (colatitude) angle at which to calculate the W function. phi : array_like Azimuthal angle at which to calculate the W function. @@ -927,9 +971,16 @@ def spin_wigner(rho, theta, phi): Values representing the spin Wigner function at the values specified by THETA and PHI. - Notes - ----- - Experimental. + References + ---------- + [1] Agarwal, G. S. (1981). Phys. Rev. A, 24(6), 2889–2896. + https://doi.org/10.1103/PhysRevA.24.2889 + + [2] Dowling, J. P., Agarwal, G. S., & Schleich, W. P. (1994). + Phys. Rev. A, 49(5), 4101–4109. https://doi.org/10.1103/PhysRevA.49.4101 + + [3] Conversion between Wigner 3-j symbol and Clebsch-Gordan coefficients + taken from Wikipedia (https://en.wikipedia.org/wiki/3-j_symbol) """ @@ -951,4 +1002,4 @@ def spin_wigner(rho, theta, phi): # sph_harm takes azimuthal angle then polar angle as arguments W += _rho_kq(rho, j, k, q) * sph_harm(q, k, PHI, THETA) - return W, THETA, PHI + return W.real, THETA, PHI