1# -*- coding: utf-8 -*- 2# This file is part of QuTiP: Quantum Toolbox in Python. 3# 4# Copyright (c) 2011 and later, Paul D. Nation and Robert J. Johansson. 5# All rights reserved. 6# 7# Redistribution and use in source and binary forms, with or without 8# modification, are permitted provided that the following conditions are 9# met: 10# 11# 1. Redistributions of source code must retain the above copyright notice, 12# this list of conditions and the following disclaimer. 13# 14# 2. Redistributions in binary form must reproduce the above copyright 15# notice, this list of conditions and the following disclaimer in the 16# documentation and/or other materials provided with the distribution. 17# 18# 3. Neither the name of the QuTiP: Quantum Toolbox in Python nor the names 19# of its contributors may be used to endorse or promote products derived 20# from this software without specific prior written permission. 21# 22# THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS 23# "AS IS" AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT 24# LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS FOR A 25# PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE COPYRIGHT 26# HOLDER OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, INCIDENTAL, 27# SPECIAL, EXEMPLARY, OR CONSEQUENTIAL DAMAGES (INCLUDING, BUT NOT 28# LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; LOSS OF USE, 29# DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED AND ON ANY 30# THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, OR TORT 31# (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT OF THE USE 32# OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF SUCH DAMAGE. 33############################################################################### 34 35import platform 36 37import numpy as np 38import pytest 39 40from qutip import ( 41 Qobj, tensor, fock_dm, basis, destroy, qdiags, sigmax, sigmay, sigmaz, 42 qeye, rand_ket, rand_super_bcsz, rand_ket_haar, rand_dm_ginibre, rand_dm, 43 rand_unitary, rand_unitary_haar, to_super, to_choi, kraus_to_choi, 44) 45from qutip.qip.operations import ( 46 hadamard_transform, swap, 47) 48# These ones are the metrics functions that we actually want to test. 49from qutip import ( 50 fidelity, tracedist, hellinger_dist, dnorm, average_gate_fidelity, 51 unitarity, hilbert_dist, bures_dist, 52) 53 54 55@pytest.fixture(scope="function", params=[2, 5, 10, 15, 25, 100]) 56def dimension(request): 57 # There are also some cases in the file where this fixture is explicitly 58 # overridden by a more local mark. That is deliberate; this dimension is 59 # intended for non-superoperators, and may cause inordinantly long tests if 60 # (for example) something uses dimension=100 then makes a superoperator out 61 # of it. 62 return request.param 63 64 65@pytest.fixture(scope="function", params=[ 66 pytest.param(rand_ket_haar, id="pure"), 67 pytest.param(rand_dm_ginibre, id="mixed"), 68]) 69def state(request, dimension): 70 return request.param(dimension) 71 72 73# Also parametrise left, right as if they're the names of two states for tests 74# that need to take two states. 75left = right = state 76 77 78# The class names have an unusual naming convention to make them more 79# convenient to use with the `pytest -k "expr"` selection syntax. They start 80# with the standard `Test`, but then are the name of the function they are 81# testing in the function naming convention, so it's easy to remember the 82# selector to choose a particular function. 83 84class Test_fidelity: 85 def test_mixed_state_inequality(self, dimension): 86 tol = 1e-7 87 rho1 = rand_dm(dimension, 0.25) 88 rho2 = rand_dm(dimension, 0.25) 89 F = fidelity(rho1, rho2) 90 assert 1 - F <= np.sqrt(1 - F*F) + tol 91 92 @pytest.mark.parametrize('right_dm', [True, False], ids=['mixed', 'pure']) 93 @pytest.mark.parametrize('left_dm', [True, False], ids=['mixed', 'pure']) 94 def test_orthogonal(self, left_dm, right_dm, dimension): 95 left = basis(dimension, 0) 96 right = basis(dimension, dimension//2) 97 if left_dm: 98 left = left.proj() 99 if right_dm: 100 right = right.proj() 101 assert fidelity(left, right) == pytest.approx(0, abs=1e-6) 102 103 def test_invariant_under_unitary_transformation(self, dimension): 104 rho1 = rand_dm(dimension, 0.25) 105 rho2 = rand_dm(dimension, 0.25) 106 U = rand_unitary(dimension, 0.25) 107 F = fidelity(rho1, rho2) 108 FU = fidelity(U*rho1*U.dag(), U*rho2*U.dag()) 109 assert F == pytest.approx(FU, rel=1e-5) 110 111 def test_state_with_itself(self, state): 112 assert fidelity(state, state) == pytest.approx(1, abs=1e-6) 113 114 def test_bounded(self, left, right, dimension): 115 """Test that fidelity is bounded on [0, 1].""" 116 tol = 1e-7 117 assert -tol <= fidelity(left, right) <= 1 + tol 118 119 def test_pure_state_equivalent_to_overlap(self, dimension): 120 """Check fidelity against pure-state overlap, see gh-361.""" 121 psi = rand_ket_haar(dimension) 122 phi = rand_ket_haar(dimension) 123 overlap = np.abs(psi.overlap(phi)) 124 assert fidelity(psi, phi) == pytest.approx(overlap, abs=1e-7) 125 126 ket_0 = basis(2, 0) 127 ket_1 = basis(2, 1) 128 ket_p = (ket_0 + ket_1).unit() 129 ket_py = (ket_0 + np.exp(0.25j*np.pi)*ket_1).unit() 130 max_mixed = qeye(2).unit() 131 132 @pytest.mark.parametrize(['left', 'right', 'expected'], [ 133 pytest.param(ket_0, ket_p, np.sqrt(0.5), id="|0>,|+>"), 134 pytest.param(ket_0, ket_1, 0, id="|0>,|1>"), 135 pytest.param(ket_0, max_mixed, np.sqrt(0.5), id="|0>,id/2"), 136 pytest.param(ket_p, ket_py, np.sqrt(0.125 + (0.5+np.sqrt(0.125))**2), 137 id="|+>,|+'>"), 138 ]) 139 def test_known_cases(self, left, right, expected): 140 assert fidelity(left, right) == pytest.approx(expected, abs=1e-7) 141 142 143class Test_tracedist: 144 def test_state_with_itself(self, state): 145 assert tracedist(state, state) == pytest.approx(0, abs=1e-6) 146 147 @pytest.mark.parametrize('right_dm', [True, False], ids=['mixed', 'pure']) 148 @pytest.mark.parametrize('left_dm', [True, False], ids=['mixed', 'pure']) 149 def test_orthogonal(self, left_dm, right_dm, dimension): 150 left = basis(dimension, 0) 151 right = basis(dimension, dimension//2) 152 if left_dm: 153 left = left.proj() 154 if right_dm: 155 right = right.proj() 156 assert tracedist(left, right) == pytest.approx(1, abs=1e-6) 157 158 def test_invariant_under_unitary_transformation(self, dimension): 159 rho1 = rand_dm(dimension, 0.25) 160 rho2 = rand_dm(dimension, 0.25) 161 U = rand_unitary(dimension, 0.25) 162 D = tracedist(rho1, rho2) 163 DU = tracedist(U*rho1*U.dag(), U*rho2*U.dag()) 164 assert D == pytest.approx(DU, rel=1e-5) 165 166 167class Test_hellinger_dist: 168 @pytest.mark.parametrize('right_dm', [True, False], ids=['mixed', 'pure']) 169 @pytest.mark.parametrize('left_dm', [True, False], ids=['mixed', 'pure']) 170 def test_orthogonal(self, left_dm, right_dm, dimension): 171 left = basis(dimension, 0) 172 right = basis(dimension, dimension//2) 173 if left_dm: 174 left = left.proj() 175 if right_dm: 176 right = right.proj() 177 expected = np.sqrt(2) 178 assert hellinger_dist(left, right) == pytest.approx(expected, abs=1e-6) 179 180 def test_state_with_itself(self, state): 181 assert hellinger_dist(state, state) == pytest.approx(0, abs=1e-6) 182 183 def test_known_cases_pure_states(self, dimension): 184 left = rand_ket(dimension) 185 right = rand_ket(dimension) 186 expected = np.sqrt(2 * (1 - np.abs(left.overlap(right))**2)) 187 assert hellinger_dist(left, right) == pytest.approx(expected, abs=1e-7) 188 189 @pytest.mark.parametrize('dimension', [2, 5, 10, 25]) 190 def test_monotonicity(self, dimension): 191 """ 192 Check monotonicity w.r.t. tensor product, see. Eq. (45) in 193 arXiv:1611.03449v2: 194 hellinger_dist(rhoA & rhoB, sigmaA & sigmaB) 195 >= hellinger_dist(rhoA, sigmaA) 196 with equality iff sigmaB = rhoB where '&' is the tensor product. 197 """ 198 tol = 1e-5 199 rhoA, rhoB, sigmaA, sigmaB = [rand_dm(dimension) for _ in [None]*4] 200 rho = tensor(rhoA, rhoB) 201 rho_sim = tensor(rhoA, sigmaB) 202 sigma = tensor(sigmaA, sigmaB) 203 dist = hellinger_dist(rhoA, sigmaA) 204 assert hellinger_dist(rho, sigma) + tol > dist 205 assert hellinger_dist(rho_sim, sigma) == pytest.approx(dist, abs=tol) 206 207 208# TODO: resolve the Mac failures. 209@pytest.mark.skipif( 210 "Darwin" in platform.system(), 211 reason="average gate fidelity tests broken on macOS as of July 2019", 212) 213class Test_average_gate_fidelity: 214 def test_identity(self, dimension): 215 id = qeye(dimension) 216 assert average_gate_fidelity(id) == pytest.approx(1, abs=1e-12) 217 218 @pytest.mark.parametrize('dimension', [2, 5, 10, 20]) 219 def test_bounded(self, dimension): 220 tol = 1e-7 221 channel = rand_super_bcsz(dimension) 222 assert -tol <= average_gate_fidelity(channel) <= 1 + tol 223 224 @pytest.mark.parametrize('dimension', [2, 5, 10, 20]) 225 def test_unitaries_equal_1(self, dimension): 226 """Tests that for random unitaries U, AGF(U, U) = 1.""" 227 tol = 1e-7 228 U = rand_unitary_haar(dimension) 229 SU = to_super(U) 230 assert average_gate_fidelity(SU, target=U) == pytest.approx(1, abs=tol) 231 232 233class Test_hilbert_dist: 234 def test_known_cases(self): 235 r1 = qdiags(np.array([0.5, 0.5, 0, 0]), 0) 236 r2 = qdiags(np.array([0, 0, 0.5, 0.5]), 0) 237 assert hilbert_dist(r1, r2) == pytest.approx(1, abs=1e-6) 238 239 240paulis = [qeye(2), sigmax(), sigmay(), sigmaz()] 241 242 243class Test_unitarity: 244 @pytest.mark.parametrize(['operator', 'expected'], [ 245 pytest.param(to_super(sigmax()), 1, id="sigmax"), 246 pytest.param(0.25 * sum(to_super(x) for x in paulis), 0, id="paulis"), 247 pytest.param(0.5 * (to_super(qeye(2)) + to_super(sigmax())), 1/3, 248 id="id+sigmax"), 249 ]) 250 def test_known_cases(self, operator, expected): 251 assert unitarity(operator) == pytest.approx(expected, abs=1e-7) 252 253 @pytest.mark.parametrize('n_qubits', [1, 2, 3, 4, 5]) 254 def test_bounded(self, n_qubits): 255 tol = 1e-7 256 operator = rand_super_bcsz(2**n_qubits) 257 assert -tol <= unitarity(operator) <= 1 + tol 258 259 260class TestComparisons: 261 """Test some known inequalities between two different metrics.""" 262 263 def test_inequality_tracedist_to_fidelity(self, left, right): 264 tol = 1e-7 265 assert 1 - fidelity(left, right) <= tracedist(left, right) + tol 266 267 def test_inequality_hellinger_dist_to_bures_dist(self, left, right): 268 tol = 1e-7 269 hellinger = hellinger_dist(left, right) 270 bures = bures_dist(left, right) 271 assert bures <= hellinger + tol 272 273 274def overrotation(x): 275 return to_super((1j * np.pi * x * sigmax() / 2).expm()) 276 277 278def had_mixture(x): 279 id_chan = to_choi(qeye(2)) 280 S_eye = to_super(id_chan) 281 S_H = to_super(hadamard_transform()) 282 return (1 - x) * S_eye + x * S_H 283 284 285def swap_map(x): 286 base = (1j * x * swap()).expm() 287 dims = [[[2], [2]], [[2], [2]]] 288 return Qobj(base, dims=dims, type='super', superrep='super') 289 290 291def adc_choi(x): 292 kraus = [ 293 np.sqrt(1 - x) * qeye(2), 294 np.sqrt(x) * destroy(2), 295 np.sqrt(x) * fock_dm(2, 0), 296 ] 297 return kraus_to_choi(kraus) 298 299 300# dnorm tests have always been slightly flaky; in some cases, cvxpy will fail 301# to solve the problem, and this can cause an entire test-suite failure. As 302# long as we are using random tests (perhaps not ideal), this will happen 303# occasionally. This isn't entirely a bug, it's just a reality of using a 304# one-size-fits-all solver; we've historically assumed users who come up 305# against this sort of thing will be accepting of the fact that dnorm 306# calculation is nontrivial, and isn't always entirely feasible. 307# 308# To deal with it, we allow each test to be rerun twice, using 309# pytest-rerunfailures. This should forbid pathological cases where the test 310# is failing every time, but not penalise one-off failures. As far as we know, 311# the failing tests always involve a random step, so triggering a re-run will 312# have them choose new variables as well. 313@pytest.mark.flaky(reruns=2) 314class Test_dnorm: 315 # Skip dnorm tests if we don't have cvxpy or cvxopt available, since it 316 # depends on them. 317 cvxpy = pytest.importorskip("cvxpy") 318 cvxopt = pytest.importorskip("cvxopt") 319 320 @pytest.fixture(params=[2, 3]) 321 def dimension(self, request): 322 return request.param 323 324 @pytest.fixture(params=[True, False], ids=['sparse', 'dense']) 325 def sparse(self, request): 326 return request.param 327 328 @pytest.mark.parametrize("variable", [0.1, 0.5, 0.9]) 329 def test_sparse_against_dense_adc(self, variable): 330 """ 331 Test sparse versus dense dnorm calculation for a sample of 332 amplitude-damping channels. 333 """ 334 # Choi matrix for identity channel on 1 qubit 335 A = kraus_to_choi([qeye(2)]) 336 B = adc_choi(variable) 337 dense = dnorm(A, B, force_solve=True, sparse=False) 338 sparse = dnorm(A, B, force_solve=True, sparse=True) 339 assert dense == pytest.approx(sparse, abs=1e-7) 340 341 @pytest.mark.repeat(3) 342 def test_sparse_against_dense_random(self, dimension): 343 """ 344 Test sparse versus dense dnorm calculation for random superoperators. 345 """ 346 A = rand_super_bcsz(dimension) 347 dense_run_result = dnorm(A, force_solve=True, sparse=False) 348 sparse_run_result = dnorm(A, force_solve=True, sparse=True) 349 assert dense_run_result == pytest.approx(sparse_run_result, abs=1e-7) 350 351 def test_bounded(self, dimension, sparse): 352 """dnorm(A - B) in [0, 2] for random superops A, B.""" 353 tol = 1e-7 354 A, B = rand_super_bcsz(dimension), rand_super_bcsz(dimension) 355 assert -tol <= dnorm(A, B, sparse=sparse) <= 2 + tol 356 357 def test_qubit_simple_known_cases(self, sparse): 358 """Check agreement for known qubit channels.""" 359 id_chan = to_choi(qeye(2)) 360 X_chan = to_choi(sigmax()) 361 depol = to_choi(Qobj( 362 qeye(4), dims=[[[2], [2]], [[2], [2]]], superrep='chi', 363 )) 364 # We need to restrict the number of iterations for things on the 365 # boundary, such as perfectly distinguishable channels. 366 assert ( 367 dnorm(id_chan, X_chan, sparse=sparse) 368 == pytest.approx(2, abs=1e-7) 369 ) 370 assert ( 371 dnorm(id_chan, depol, sparse=sparse) 372 == pytest.approx(1.5, abs=1e-7) 373 ) 374 # Finally, we add a known case from Johnston's QETLAB documentation, 375 # || Phi - I ||_♢, 376 # where Phi(X) = UXU⁺ and U = [[1, 1], [-1, 1]] / sqrt(2). 377 U = np.sqrt(0.5) * Qobj([[1, 1], [-1, 1]]) 378 assert ( 379 dnorm(U, qeye(2), sparse=sparse) 380 == pytest.approx(np.sqrt(2), abs=1e-7) 381 ) 382 383 @pytest.mark.parametrize(["variable", "expected", "generator"], [ 384 [1.0e-3, 3.141591e-03, overrotation], 385 [3.1e-3, 9.738899e-03, overrotation], 386 [1.0e-2, 3.141463e-02, overrotation], 387 [3.1e-2, 9.735089e-02, overrotation], 388 [1.0e-1, 3.128689e-01, overrotation], 389 [3.1e-1, 9.358596e-01, overrotation], 390 [1.0e-3, 2.000000e-03, had_mixture], 391 [3.1e-3, 6.200000e-03, had_mixture], 392 [1.0e-2, 2.000000e-02, had_mixture], 393 [3.1e-2, 6.200000e-02, had_mixture], 394 [1.0e-1, 2.000000e-01, had_mixture], 395 [3.1e-1, 6.200000e-01, had_mixture], 396 [1.0e-3, 2.000000e-03, swap_map], 397 [3.1e-3, 6.199997e-03, swap_map], 398 [1.0e-2, 1.999992e-02, swap_map], 399 [3.1e-2, 6.199752e-02, swap_map], 400 [1.0e-1, 1.999162e-01, swap_map], 401 [3.1e-1, 6.173918e-01, swap_map], 402 ]) 403 def test_qubit_known_cases(self, variable, expected, generator, sparse): 404 """ 405 Test cases based on comparisons to pre-existing dnorm implementations. 406 In particular, the targets for the following test cases were generated 407 using QuantumUtils for MATLAB (https://goo.gl/oWXhO9). 408 """ 409 id_chan = to_choi(qeye(2)) 410 channel = generator(variable) 411 assert ( 412 dnorm(channel, id_chan, sparse=sparse) 413 == pytest.approx(expected, abs=1e-7) 414 ) 415 416 def test_qubit_scalar(self, dimension): 417 """dnorm(a * A) == a * dnorm(A) for scalar a, qobj A.""" 418 a = np.random.random() 419 A = rand_super_bcsz(dimension) 420 B = rand_super_bcsz(dimension) 421 assert dnorm(a*A, a*B) == pytest.approx(a*dnorm(A, B), abs=1e-7) 422 423 def test_qubit_triangle(self, dimension): 424 """Check that dnorm(A + B) <= dnorm(A) + dnorm(B).""" 425 A = rand_super_bcsz(dimension) 426 B = rand_super_bcsz(dimension) 427 assert dnorm(A + B) <= dnorm(A) + dnorm(B) + 1e-7 428 429 @pytest.mark.repeat(3) 430 @pytest.mark.parametrize("generator", [ 431 pytest.param(rand_super_bcsz, id="super"), 432 pytest.param(rand_unitary_haar, id="unitary"), 433 ]) 434 def test_force_solve(self, dimension, generator): 435 """ 436 Metrics: checks that special cases for dnorm agree with SDP solutions. 437 """ 438 A, B = generator(dimension), generator(dimension) 439 assert ( 440 dnorm(A, B, force_solve=False) 441 == pytest.approx(dnorm(A, B, force_solve=True), abs=1e-5) 442 ) 443 444 @pytest.mark.repeat(3) 445 def test_cptp(self, dimension, sparse): 446 """Check that the diamond norm is one for CPTP maps.""" 447 A = rand_super_bcsz(dimension) 448 assert A.iscptp 449 assert dnorm(A, sparse=sparse) == pytest.approx(1, abs=1e-7) 450