# encoding: utf-8
from __future__ import division, print_function
import itertools as it
from inspect import isfunction
import mpnum as mp
import mpnum.factory as factory
import mpnum.mpsmpo as mpsmpo
import mpnum.povm as povm
import numpy as np
import pytest as pt
from _pytest.mark import matchmark
from mpnum import utils
from mpnum.utils.pmf import project_pmf
from numpy.testing import assert_almost_equal, assert_array_almost_equal
from six.moves import range, zip, zip_longest
ALL_POVMS = {name: constructor for name, constructor in povm.__dict__.items()
if name.endswith('_povm') and isfunction(constructor)}
# nr_sites, startsite, local_dim
MPPOVM_PARAM = [
(4, 0, 2), (5, 0, 3), (5, 1, 2), (6, 1, 2),
pt.mark.long((7, 1, 2)), pt.mark.long((8, 2, 2)), pt.mark.long((6, 1, 3))
]
# --------------------------------------------------------------
# Many tests in this module have a NON-ZERO FAILURE PROBABILITY.
# --------------------------------------------------------------
#
# It is not guaranteed that probabilities will be estimated to within
# the thresholds used in the tests. Exceptions should have a low
# probability.
#
# TODO: Many tests in this module basically check that "estimation
# error becomes smaller as the number of samples increases". Usually,
# we check that the estimation error is smaller than a constant times
# 1/sqrt(number of samples) (cf. e.g. central limit theorem or
# variance of the estimator). This is a basic consistency check but
# far from a statistical correctness check, which would be desirable.
#
# method, n_samples
MPPOVM_SAMPLE_PARAM = [
('direct', 100), ('cond', 100), pt.mark.long(('cond', 1000)),
('direct', 2500), pt.mark.long(('direct', 10000)),
pt.mark.long(('direct', 40000)), pt.mark.long(('direct', 80000))
]
def mp_from_array_repeat(array, nr_sites):
"""Generate a MPA representation of the `nr_sites`-fold tensor product of
array.
"""
mpa = mp.MPArray.from_array(array)
return mp.chain(it.repeat(mpa, nr_sites))
@pt.fixture(params=['random', 'pauli'])
def nopovm(request, local_dim, rgen):
"""Provide different POVMs and non-POVMs for testing
We provide instances of :class:`povm.localpovm.POVM` with the
following elements:
* `pauli`: Generated by :func:`povm.pauli_povm()`
* `random`: Random (non-Hermitian, non-positive) elements for
testing. (These elements do not constitute a POVM. We use them
to distinguish elem.conj() from elem.T in our code.)
"""
nopovm_name = request.param
if nopovm_name == 'pauli':
return povm.pauli_povm(local_dim)
elif nopovm_name == 'random':
d = local_dim
return povm.localpovm.POVM(factory._zrandn((2 * d**2, d, d), rgen))
else:
raise ValueError('Unknown fixture name {}'.format(nopovm_name))
@pt.mark.parametrize('dim', [(2), (3), (6), (7)])
def test_povm_normalization_ic(dim):
for name, constructor in ALL_POVMS.items():
# Check that the POVM is normalized: elements must sum to the identity
current_povm = constructor(dim)
element_sum = sum(iter(current_povm))
assert_array_almost_equal(element_sum, np.eye(dim))
# Check that the attribute that says whether the POVM is IC is correct.
linear_inversion_recons = np.dot(current_povm.linear_inversion_map,
current_povm.probability_map)
if current_povm.informationally_complete:
assert_array_almost_equal(
linear_inversion_recons, np.eye(dim**2),
err_msg='POVM {} is not informationally complete'.format(name))
else:
assert np.abs(linear_inversion_recons - np.eye(dim**2)).max() > 0.1, \
'POVM {} is informationally complete'.format(name)
@pt.mark.parametrize('nr_sites, local_dim, rank',
[(6, 2, 7), (3, 3, 3), (3, 6, 3), (3, 7, 4)])
def test_povm_ic_mpa(nr_sites, local_dim, rank, rgen):
# Check that the tensor product of the PauliGen POVM is IC.
paulis = povm.pauli_povm(local_dim)
inv_map = mp_from_array_repeat(paulis.linear_inversion_map, nr_sites)
probab_map = mp_from_array_repeat(paulis.probability_map, nr_sites)
reconstruction_map = mp.dot(inv_map, probab_map)
eye = factory.eye(nr_sites, local_dim**2)
assert mp.norm(reconstruction_map - eye) < 1e-5
# Check linear inversion for a particular example MPA.
# Linear inversion works for arbitrary matrices, not only for states,
# so we test it for an arbitrary MPA.
# Normalize, otherwise the absolute error check below will not work.
mpa = factory.random_mpa(nr_sites, local_dim**2, rank,
dtype=np.complex_, randstate=rgen, normalized=True)
probabs = mp.dot(probab_map, mpa)
recons = mp.dot(inv_map, probabs)
assert mp.norm(recons - mpa) < 1e-6
@pt.mark.parametrize('local_dim', [(2), (3), (6), (7)])
def test_povm_probability_map(local_dim, nopovm, rgen):
# Use a random matrix rho for testing (instead of a positive matrix).
rho = factory._zrandn((local_dim, local_dim), rgen)
# Compare output from `povm.localpovm.POVM.probability_map` with
# calculating probabilities element by element.
probab_direct = np.array([np.trace(np.dot(elem, rho)) for elem in nopovm])
probab_pmap = np.dot(nopovm.probability_map, rho.ravel())
assert_array_almost_equal(probab_pmap, probab_direct)
@pt.mark.parametrize('nr_sites, width, local_dim, rank',
[(6, 3, 2, 5), (4, 2, 3, 4)])
def test_mppovm_expectation(nr_sites, width, local_dim, rank, nopovm, rgen):
# Verify that :func:`povm.MPPovm.expectations()` produces
# correct results.
pmap = nopovm.probability_map
mpnopovm = povm.MPPovm.from_local_povm(nopovm, width)
# Use a random MPO rho for testing (instead of a positive MPO).
rho = factory.random_mpo(nr_sites, local_dim, rank, rgen)
reductions = mpsmpo.reductions_mpo(rho, width)
# Compute expectation values with mpnopovm.expectations(), which
# uses mpnopovm.probability_map.
expectations = list(mpnopovm.expectations(rho))
assert len(expectations) == nr_sites - width + 1
for evals_mp, rho_red in zip_longest(expectations, reductions):
# Compute expectation values by constructing each tensor
# product POVM element.
rho_red_matrix = rho_red.to_array_global().reshape(
(local_dim**width,) * 2)
evals = []
for factors in it.product(nopovm, repeat=width):
elem = utils.mkron(*factors)
evals.append(np.trace(np.dot(elem, rho_red_matrix)))
evals = np.array(evals).reshape((len(nopovm),) * width)
# Compute expectation with a different construction. In the
# end, this is (should be, we verify it here) equivalent to
# what `mpnopovm.expectations()` does.
evals_ten = rho_red.ravel().to_array()
for _ in range(width):
evals_ten = np.tensordot(evals_ten, pmap, axes=(0, 1))
assert_array_almost_equal(evals_ten, evals)
assert_array_almost_equal(evals_mp.to_array(), evals)
@pt.mark.parametrize(
'nr_sites, local_dim, rank, startsite, width',
[(4, 2, 3, 0, 4), (7, 2, 3, 1, 3), (6, (7, 3, 2, 5, 2, 3), 3, 2, 3)]
)
def test_mppovm_embed_expectation(
nr_sites, local_dim, rank, startsite, width, rgen):
if hasattr(local_dim, '__iter__'):
local_dim2 = local_dim
else:
local_dim2 = [local_dim] * nr_sites
local_dim2 = list(zip(local_dim2, local_dim2))
# Create a local POVM `red_povm`, embed it onto a larger chain
# (`full_povm`), and go back to the reduced POVM.
red_povm = mp.chain(
mp.povm.MPPovm.from_local_povm(mp.povm.pauli_povm(d), 1)
for d, _ in local_dim2[startsite:startsite + width]
)
full_povm = red_povm.embed(nr_sites, startsite, local_dim)
axes = [(1, 2) if i < startsite or i >= startsite + width else None
for i in range(nr_sites)]
red_povm2 = mp.partialtrace(full_povm, axes, mp.MPArray)
red_povm2 = mp.prune(red_povm2, singletons=True)
red_povm2 /= np.prod([d for i, (d, _) in enumerate(local_dim2)
if i < startsite or i >= startsite + width])
assert_almost_equal(mp.normdist(red_povm, red_povm2), 0.0)
# Test with an arbitrary random MPO instead of an MPDO
mpo = mp.factory.random_mpa(nr_sites, local_dim2, rank, rgen,
dtype=np.complex_, normalized=True)
mpo_red = next(mp.reductions_mpo(mpo, width, startsites=[startsite]))
ept = mp.prune(full_povm.pmf(mpo, 'mpdo'), singletons=True).to_array()
ept_red = red_povm.pmf(mpo_red, 'mpdo').to_array()
assert_array_almost_equal(ept, ept_red)
@pt.mark.parametrize('nr_sites, width, local_dim, rank',
[(6, 3, 2, 5), (4, 2, 3, 4)])
def test_mppovm_expectation_pure(nr_sites, width, local_dim, rank, rgen):
paulis = povm.pauli_povm(local_dim)
mppaulis = povm.MPPovm.from_local_povm(paulis, width)
psi = factory.random_mps(nr_sites, local_dim, rank, randstate=rgen)
rho = mpsmpo.mps_to_mpo(psi)
expect_psi = list(mppaulis.expectations(psi))
expect_rho = list(mppaulis.expectations(rho))
assert len(expect_psi) == len(expect_rho)
for e_rho, e_psi in zip(expect_rho, expect_psi):
assert_array_almost_equal(e_rho.to_array(), e_psi.to_array())
@pt.mark.parametrize('nr_sites, width, local_dim, rank',
[(6, 3, 2, 5), (4, 2, 3, 4)])
def test_mppovm_expectation_pmps(nr_sites, width, local_dim, rank, rgen):
paulis = povm.pauli_povm(local_dim)
mppaulis = povm.MPPovm.from_local_povm(paulis, width)
pmps = factory.random_mpa(nr_sites, (local_dim, local_dim), rank,
dtype=np.complex_, randstate=rgen)
rho = mpsmpo.pmps_to_mpo(pmps)
expect_psi = list(mppaulis.expectations(pmps, mode='pmps'))
expect_rho = list(mppaulis.expectations(rho))
assert len(expect_psi) == len(expect_rho)
for e_rho, e_psi in zip(expect_rho, expect_psi):
assert_array_almost_equal(e_rho.to_array(), e_psi.to_array())
@pt.mark.parametrize(
'nr_sites, local_dim, rank, startsite, width',
[(4, 2, 3, 0, 4), (4, ((5, 2), (2, 3), (3, 2), (2, 2)), 3, 0, 4),
(7, 2, 3, 1, 3), (6, ((5, 2), (2, 3), (3, 2), (2, 2), (5, 3), (3, 2)), 3, 2, 3)]
)
def test_mppovm_pmf_as_array_pmps(
nr_sites, local_dim, rank, startsite, width, rgen):
if hasattr(local_dim, '__len__'):
pdims = [d for d, _ in local_dim]
mppaulis = mp.chain(
povm.MPPovm.from_local_povm(povm.pauli_povm(d), 1)
for d in pdims[startsite:startsite + width]
)
else:
pdims = local_dim
local_dim = (local_dim, local_dim)
mppaulis = povm.MPPovm.from_local_povm(povm.pauli_povm(pdims), width)
mppaulis = mppaulis.embed(nr_sites, startsite, pdims)
pmps = factory.random_mpa(nr_sites, local_dim, rank,
dtype=np.complex_, randstate=rgen, normalized=True)
rho = mpsmpo.pmps_to_mpo(pmps)
expect_rho = mppaulis.pmf_as_array(rho, 'mpdo')
for impl in ['default', 'pmps-ltr', 'pmps-symm']:
expect_pmps = mppaulis.pmf_as_array(pmps, 'pmps', impl=impl)
assert_array_almost_equal(expect_rho, expect_pmps, err_msg=impl)
@pt.mark.benchmark(group='pmf_as_array_pmps')
@pt.mark.parametrize(
'nr_sites, local_dim, rank, startsite, width', [(10, 2, 16, 0, 10)])
@pt.mark.parametrize('impl', ['default', 'pmps-ltr', 'pmps-symm'])
def test_mppovm_pmf_as_array_pmps_benchmark(
nr_sites, local_dim, rank, startsite, width, impl, rgen, benchmark):
pauli_y = povm.pauli_parts(local_dim)[1]
mpp_y = povm.MPPovm.from_local_povm(pauli_y, width) \
.embed(nr_sites, startsite, local_dim)
pmps = factory.random_mpa(nr_sites, (local_dim, local_dim), rank,
dtype=np.complex_, randstate=rgen, normalized=True)
benchmark(lambda: mpp_y.pmf_as_array(pmps, 'pmps', impl=impl))
@pt.mark.benchmark(group='pmf_as_array_pmps2', min_rounds=2)
@pt.mark.parametrize(
'nr_sites, local_dim, rank, startsite, width',
[(32, 2, 20, 10, 6)])
@pt.mark.parametrize('impl', ['default', 'pmps-ltr', 'pmps-symm'])
def test_mppovm_pmf_as_array_pmps_benchmark2(
nr_sites, local_dim, rank, startsite, width, impl, rgen, benchmark):
return test_mppovm_pmf_as_array_pmps_benchmark(
nr_sites, local_dim, rank, startsite, width, impl, rgen, benchmark)
@pt.mark.parametrize(
'nr_sites, n_small, small_startsite, local_dim',
[(5, 2, 1, 2), (5, 2, 0, 2), (5, 2, 3, 2),
(8, 3, 2, 2), (8, 3, 2, 3), (40, 3, 10, 2)])
def test_mppovm_match_elems_local(
nr_sites, n_small, small_startsite, local_dim, eps=1e-10):
"""Check that match_elems() works for single- and
multi-Pauli MPPOVMs"""
n_right = nr_sites - n_small - small_startsite
assert n_right >= 0
# "Big" POVM: X on all sites, "small" POVM: X on `n_small` neighbours
x = povm.x_povm(local_dim)
big = povm.MPPovm.from_local_povm(x, nr_sites)
small = povm.MPPovm.from_local_povm(x, n_small) \
.embed(nr_sites, small_startsite, local_dim)
match, prefactors = small.match_elems(big)
assert match.shape == tuple([len(x)] * n_small * 2)
assert match.shape == prefactors.shape
# Verify the expected one-to-one correspondence between POVM elements.
want = np.eye(np.prod(small.outdims), dtype=bool).reshape(match.shape)
assert (match == want).all()
assert (abs(prefactors[match] - 1.0) <= eps).all()
assert np.isnan(prefactors[~match]).all()
# "Big" POVM: X on all sites, "small" POVM: Paulis on `n_small` neighbours
paulis = povm.pauli_povm(local_dim)
small = povm.MPPovm.from_local_povm(paulis, n_small) \
.embed(nr_sites, small_startsite, local_dim)
match, prefactors = small.match_elems(big)
assert match.shape == tuple([len(paulis)] * n_small + [len(x)] * n_small)
assert match.shape == prefactors.shape
# Verify that the X POVM elements were found where we expect them
x_pos = tuple([slice(0, len(x))] * n_small)
want = np.zeros_like(match)
want[x_pos] = np.eye(len(x)**n_small, dtype=bool).reshape([len(x)] * 2 * n_small)
assert (match == want).all()
want = (2 if local_dim > 2 else 3)**-n_small
assert (abs(prefactors[match] - want) / want <= eps).all()
assert np.isnan(prefactors[~match]).all()
# "Big" POVM: Y on all sites, "small" POVM: Paulis on `n_small` neighbours
y = povm.y_povm(local_dim)
big = povm.MPPovm.from_local_povm(y, nr_sites)
match, prefactors = small.match_elems(big)
assert match.shape == tuple([len(paulis)] * n_small + [len(y)] * n_small)
assert match.shape == prefactors.shape
# Verify that the Y POVM elements were found where we expect them
y_pos = tuple([slice(len(x), len(x) + len(y))] * n_small)
want = np.zeros_like(match)
want[y_pos] = np.eye(len(y)**n_small, dtype=bool).reshape([len(y)] * 2 * n_small)
assert (match == want).all()
want = (2 if local_dim > 2 else 3)**-n_small
assert (abs(prefactors[match] - want) / want <= eps).all()
assert np.isnan(prefactors[~match]).all()
def test_mppovm_match_elems_bell(eps=1e-10):
"""Test match_elems() for a non-product MPPovm"""
# Four Bell states (basis: |00>, |01>, |10>, |11>)
bell = np.array((
[(1/3)**0.5, 0, 0, (1/3)**0.5], # (0, 0): |00> + |11> (proj. weight 1/3)
[0, 1, 1, 0], # (0, 1): |01> + |10> (proj. weight 1)
[(2/3)**0.5, 0, 0, -(2/3)**0.5], # (0, 2): |00> - |11> (proj. weight 2/3)
[0, 1, -1, 0], # (1, 0): |01> - |10> (proj. weight 1)
[(1/3)**0.5, 0, 0, -(1/3)**0.5], # (1, 1): |00> - |11> (proj. weight 1/3)
[(2/3)**0.5, 0, 0, (2/3)**0.5], # (1, 2): |00> + |11> (proj. weight 2/3)
)) / 2**0.5
bell_proj = np.einsum('ij, ik -> ijk', bell, bell.conj())
bell_proj = bell_proj.reshape((2, 3) + (2,) * 4)
# Four Bell states and two product states
vecs = np.array((
[0, 1, -1, 0], # (0, 0): |01> - |10> (proj. weight 0.5)
[0, 2**0.5, 0, 0], # (0, 1): |01> (proj. weight 0.5)
[0, 0, 2**0.5, 0], # (1, 0): |10> (proj. weight 0.5)
[2**0.5, 0, 0, -2**0.5], # (1, 1): |00> - |11> (proj. weight 1)
[0, 1, 1, 0], # (2, 0): |01> + |10> (proj. weight 0.5)
[2**0.5, 0, 0, 2**0.5], # (2, 1): |00> + |11> (proj. weight 1)
)) / 2
proj = np.einsum('ij, ik -> ijk', vecs, vecs.conj())
proj = proj.reshape((3,) + (2,) * 5)
# Big POVM: The four Bell states (repeated two times)
big = povm.MPPovm.from_array_global(bell_proj, ndims=3)
big = mp.chain([big, big])
# Small POVM: Two of the Bell states and four product states (on
# the last two sites)
small = povm.MPPovm.from_array_global(proj, ndims=3).embed(4, 2, 2)
# Check that the POVM is normalized: elements must sum to the identity
for mppovm in big, small:
element_sum = sum(x.to_array_global().reshape(16, 16)
for x in mppovm.elements)
assert_array_almost_equal(element_sum, np.eye(16))
match, prefactors = small.match_elems(big, eps=eps)
# Verify the correspondence which can be read off above
want = np.zeros((3, 2, 2, 3), dtype=bool)
want[0, 0, 1, 0] = True # |01> - |10>
want[2, 0, 0, 1] = True # |01> + |10>
want[1, 1, 0, 2] = True # |00> - |11>
want[1, 1, 1, 1] = True # |00> - |11>
want[2, 1, 0, 0] = True # |00> + |11>
want[2, 1, 1, 2] = True # |00> + |11>
assert (match == want).all()
assert abs(prefactors[0, 0, 1, 0] - 0.5) <= eps
assert abs(prefactors[2, 0, 0, 1] - 0.5) <= eps
assert abs(prefactors[1, 1, 0, 2] - 1.5) <= eps
assert abs(prefactors[1, 1, 1, 1] - 3) <= eps
assert abs(prefactors[2, 1, 0, 0] - 3) <= eps
assert abs(prefactors[2, 1, 1, 2] - 1.5) <= eps
assert np.isnan(prefactors[~match]).all()
@pt.mark.parametrize('method, n_samples',
MPPOVM_SAMPLE_PARAM + [pt.mark.verylong(('cond', 10000))])
@pt.mark.parametrize('nr_sites, startsite, local_dim', MPPOVM_PARAM)
def test_mppovm_sample(
method, n_samples, nr_sites, startsite, local_dim, rgen):
"""Check that probability estimates from samples are reasonable accurate"""
rank = 3
eps = 1e-10
mps = factory.random_mps(nr_sites, local_dim, rank, rgen)
mps.canonicalize()
local_x = povm.x_povm(local_dim)
local_y = povm.y_povm(local_dim)
xx = povm.MPPovm.from_local_povm(local_x, 2)
y = povm.MPPovm.from_local_povm(local_y, 1)
mpp = mp.chain([xx, povm.MPPovm.eye([local_dim]), y]) \
.embed(nr_sites, startsite, local_dim)
pmf_exact = mpp.pmf_as_array(mps, 'mps', eps)
if n_samples > 100:
n_gr = 5
elif local_dim == 3:
n_gr = 2
else:
n_gr = 3
samples = mpp.sample(rgen, mps, n_samples, method, n_gr, 'mps', eps=eps)
pmf_est = mpp.est_pmf(samples)
assert abs(pmf_est.sum() - 1.0) <= eps
assert abs(pmf_exact - pmf_est).max() <= 3 / n_samples**0.5
@pt.mark.parametrize('method, n_samples', MPPOVM_SAMPLE_PARAM)
@pt.mark.parametrize('nr_sites, startsite, local_dim', MPPOVM_PARAM)
def test_mppovm_est_pmf_from(
method, n_samples, nr_sites, startsite, local_dim, rgen):
"""Check that probability estimates from samples are reasonable accurate"""
rank = 3
eps = 1e-10
mps = factory.random_mps(nr_sites, local_dim, rank, rgen)
mps.canonicalize()
lx = povm.x_povm(local_dim)
ly = povm.y_povm(local_dim)
lp = povm.pauli_povm(local_dim)
x = povm.MPPovm.from_local_povm(lx, 1)
y = povm.MPPovm.from_local_povm(ly, 1)
pauli = povm.MPPovm.from_local_povm(lp, 1)
xy = mp.chain((x, y))
mpp = mp.chain((xy,) * (nr_sites // 2))
if (nr_sites % 2) == 1:
mpp = mp.chain((mpp, x))
small_mpp = mp.chain((pauli, povm.MPPovm.eye([local_dim]), pauli, pauli)) \
.embed(nr_sites, startsite, local_dim)
x_given = np.arange(len(lp)) < len(lx)
y_given = ((np.arange(len(lp)) >= len(lx))
& (np.arange(len(lp)) < len(lx) + len(ly)))
given_sites = [x_given if ((startsite + i) % 2) == 0 else y_given
for i in (0, 2, 3)]
given_expected = np.einsum('i, j, k -> ijk', *given_sites)
pmf_exact = small_mpp.pmf_as_array(mps, 'mps', eps)
if n_samples > 100:
n_gr = 5
elif local_dim == 3:
n_gr = 2
else:
n_gr = 3
samples = mpp.sample(rgen, mps, n_samples, method, n_gr, 'mps', eps=eps)
est_pmf, est_n_samples = small_mpp.est_pmf_from(mpp, samples)
# In this case, we use all the samples from `mpp`.
assert est_n_samples == n_samples
given = ~np.isnan(est_pmf)
assert (given == given_expected).all()
assert abs(pmf_exact[given].sum() - est_pmf[given].sum()) <= eps
assert abs(pmf_exact[given] - est_pmf[given]).max() <= 1 / n_samples**0.5
@pt.mark.parametrize('method, n_samples', MPPOVM_SAMPLE_PARAM)
@pt.mark.parametrize('nr_sites, startsite, local_dim', MPPOVM_PARAM)
def test_mppovm_est(
method, n_samples, nr_sites, startsite, local_dim, rgen):
"""Check that estimates from .est_pmf() and .est_lfun() are reasonably
accurate
"""
rank = 3
eps = 1e-10
mps = factory.random_mps(nr_sites, local_dim, rank, rgen)
mps.canonicalize()
local_x = povm.x_povm(local_dim)
local_y = povm.y_povm(local_dim)
xx = povm.MPPovm.from_local_povm(local_x, 2)
y = povm.MPPovm.from_local_povm(local_y, 1)
mpp = mp.chain([xx, povm.MPPovm.eye([local_dim]), y]) \
.embed(nr_sites, startsite, local_dim)
p_exact = mpp.pmf_as_array(mps, 'mps', eps)
p_exact = project_pmf(p_exact, eps, eps)
cov_p_exact = np.diag(p_exact.flat) - np.outer(p_exact.flat, p_exact.flat)
samples = mpp.sample(rgen, mps, n_samples, method, 4, 'mps', eps=eps)
p_est = mpp.est_pmf(samples)
ept, cov = mpp.est_lfun(None, None, samples, None, eps)
ept_ex, single_cov_ex = mpp.lfun(None, None, mps, 'mps', eps)
# The two exact values must match
assert abs(ept_ex - p_exact.ravel()).max() <= eps
# The two exact values must match
assert abs(cov_p_exact - single_cov_ex).max() <= eps
# The two estimates must match. This verifies that we have chosen
# our estimator will be unbiased. (There are many other things we
# might want to know about our estimator.)
assert (ept == p_est.ravel()).all()
# The estimate must be close to the true value
assert abs(p_exact - p_est).max() <= 3 / n_samples**0.5
cov_ex = cov_p_exact / n_samples
# The covariances of the sample means (which we estimate here)
# decrease by 1/n_samples, so we multiply with n_samples before
# comparing to the rule-of-thumb for the estimation error.
assert abs(cov - cov_ex).max() * n_samples <= 1 / n_samples**0.5
funs = []
nsoutdims = mpp.nsoutdims
out = np.unravel_index(range(np.prod(nsoutdims)), nsoutdims)
out = np.array(out).T[:, None, :].copy()
for ind in range(np.prod(nsoutdims)):
funs.append(lambda s, ind=ind: (s == out[ind]).all(1))
# All probabilities sum to one, and we can estimate that well.
coeff = np.ones(len(funs), dtype=float)
# Test with dummy weights
weights = np.ones(n_samples, dtype=float)
sum_ept, sum_var = mpp.est_lfun(coeff, funs, samples, weights, eps)
assert abs(sum_ept - 1.0) <= eps
assert sum_var <= eps
# Check a sum of probabilities with varying signs.
coeff = ((-1)**rgen.choice(2, len(funs))).astype(float)
sum_ept, sum_var = mpp.est_lfun(coeff, funs, samples, None, eps)
ex_sum = np.inner(coeff, p_exact.flat)
ex_var = np.inner(coeff, np.dot(cov_ex, coeff))
assert abs(sum_ept - ex_sum) <= 5 / n_samples**0.5
assert abs(sum_var - ex_var) * n_samples <= 5 / n_samples**0.5
# Convert samples to counts and test again
counts = mpp.est_pmf(samples, normalize=False, eps=eps)
assert counts.sum() == n_samples
count_samples = np.array(np.unravel_index(range(np.prod(mpp.nsoutdims)),
mpp.nsoutdims)).T
weights = counts.ravel()
sum_ept2, sum_var2 = mpp.est_lfun(coeff, funs, count_samples, weights, eps)
assert abs(sum_ept - sum_ept2) <= eps
assert abs(sum_var - sum_var2) <= eps
@pt.mark.parametrize(
'method, n_samples', [
('direct', 100),
])
@pt.mark.parametrize(
'nr_sites, local_dim, rank, measure_width', [
(3, 3, 2, 2),
(3, 2, 3, 3),
(5, 2, 3, 2),
])
def test_mppovmlist_pack_unpack_samples(
method, n_samples, nr_sites, local_dim, rank, measure_width,
rgen, eps=1e-10):
"""Check that packing and unpacking samples does not change them"""
mps = factory.random_mps(nr_sites, local_dim, rank, rgen)
mps.canonicalize()
s_povm = povm.pauli_mpp(measure_width, local_dim).block(nr_sites)
samples = tuple(s_povm.sample(
rgen, mps, n_samples, method, mode='mps', pack=False, eps=eps))
packed = tuple(s_povm.pack_samples(samples))
unpacked = tuple(s_povm.unpack_samples(packed))
assert all(s.dtype == np.uint8 for s in samples)
assert all(s.dtype == np.uint8 for s in unpacked)
assert all((s == u).all() for s, u in zip(samples, unpacked))
def _pytest_want_verylong(request):
# Is there a better way to find out whether items marked with
# `long` should be run or not?
class dummy:
keywords = {'verylong': pt.mark.verylong}
return matchmark(dummy, request.config.option.markexpr)
@pt.fixture(params=[False, True])
def splitpauli(n_samples, nonuniform, request):
# We use this fixture to skip certain value combinations for
# non-long tests.
#
# Is there a better way to select certain value combinations from
# the different pt.mark.parametrize() decorators except for
# writing down all combinations by hand?
splitpauli = request.param
if (not splitpauli) or _pytest_want_verylong(request) \
or (n_samples >= 10000 and nonuniform):
return splitpauli
pt.skip("Should only be run in long tests")
return None
@pt.mark.parametrize(
'method, n_samples', [
('direct', 1000), ('direct', 10000),
pt.mark.verylong(('direct', 100000)), pt.mark.verylong(('cond', 100))
])
@pt.mark.parametrize(
'nr_sites, local_dim, rank, measure_width, local_width', [
(4, 2, 3, 2, 2),
pt.mark.verylong((5, 2, 2, 3, 2)),
pt.mark.verylong((4, 3, 2, 2, 2)),
])
@pt.mark.parametrize('nonuniform', [False, True])
def test_mppovmlist_est_pmf_from(
method, n_samples, nr_sites, local_dim, rank, measure_width,
local_width, nonuniform, splitpauli, rgen, eps=1e-10):
"""Verify that estimated probabilities from MPPovmList.est_pmf_from()
are reasonable accurate
"""
mps = factory.random_mps(nr_sites, local_dim, rank, rgen)
mps.canonicalize()
x, y = (povm.MPPovm.from_local_povm(p, 1)
for p in povm.pauli_parts(local_dim)[:2])
# POVM list with global support
g_povm = povm.pauli_mpps(measure_width, local_dim).repeat(nr_sites)
if nonuniform:
add_povm = mp.chain((nr_sites - 1) * (x,) + (y,))
g_povm = povm.MPPovmList(g_povm.mpps + (add_povm,))
# POVM list with local support
l_povm = povm.pauli_mpps if splitpauli else povm.pauli_mpp
l_povm = l_povm(local_width, local_dim).block(nr_sites)
samples = tuple(g_povm.sample(
rgen, mps, n_samples, method, mode='mps', eps=eps))
est_prob, n_samples = zip(*l_povm.est_pmf_from(g_povm, samples, eps))
exact_prob = tuple(l_povm.pmf_as_array(mps, 'mps', eps))
# Consistency check on n_samples: All entries should be equal
# unless `nonuniform` is True.
all_n_sam = np.concatenate(n_samples)
assert (not (all_n_sam == all_n_sam[0]).all()) == nonuniform
for n_sam, est, exact, mpp in zip(
n_samples, est_prob, exact_prob, l_povm.mpps):
assert est.shape == mpp.nsoutdims
assert est.shape == exact.shape
assert n_sam.shape == exact.shape
# Compare against exact probabilities
assert (abs(est - exact) / (3 / n_sam**0.5)).max() <= 1
def _get_povm(name, nr_sites, local_dim, local_width):
if name == 'global':
return povm.pauli_mpps(local_width, local_dim).repeat(nr_sites)
elif name == 'splitpauli':
return povm.pauli_mpps(local_width, local_dim).block(nr_sites)
elif name == 'pauli':
return povm.pauli_mpp(local_width, local_dim).block(nr_sites)
elif name == "all-y":
return povm.MPPovmList([povm.MPPovm.from_local_povm(
povm.pauli_parts(local_dim)[1], nr_sites)])
elif name == "local-x":
return povm.MPPovmList([
povm.MPPovm.from_local_povm(
povm.pauli_parts(local_dim)[0], local_width)
.embed(nr_sites, 0, local_dim)
])
else:
raise ValueError('Unknown MP-POVM list {!r}'.format(name))
POVM_COMBOS = [
('global', 'pauli'), ('splitpauli', 'pauli'), ('pauli', 'pauli'),
('global', 'all-y'), ('all-y', 'local-x'), ('all-y', 'pauli'),
pt.mark.verylong(('splitpauli', 'splitpauli')), ('pauli', 'splitpauli')
]
POVM_IDS = ['+'.join(getattr(x, 'args', (x,))[0]) for x in POVM_COMBOS]
@pt.fixture(params=POVM_COMBOS, ids=POVM_IDS)
def povm_combo(function, request):
# We use this fixture to skip certain value combinations for
# non-long tests.
#
# Is there a better way to select certain value combinations from
# the different pt.mark.parametrize() decorators except for
# writing down all combinatiosn by hand?
combo = request.param
if _pytest_want_verylong(request):
return combo
if function == 'randn' or combo == ('global', 'pauli'):
return combo
pt.skip("Should only be run in long tests")
return None
@pt.mark.parametrize(
'method, n_samples', [
pt.mark.verylong(('cond', 100)),
('direct', 1000),
pt.mark.verylong(('direct', 100000)),
])
@pt.mark.parametrize(
'nr_sites, local_dim, rank, measure_width, local_width', [
(3, 2, 3, 2, 2),
pt.mark.verylong((4, 2, 3, 3, 2)),
pt.mark.verylong((5, 2, 3, 2, 2)),
])
@pt.mark.parametrize('nonuniform', [True, pt.mark.verylong(False)])
@pt.mark.parametrize('function',
['randn', 'ones', 'signs', pt.mark.verylong('rand')])
def test_mppovmlist_est_lfun_from(
method, n_samples, nr_sites, local_dim, rank, measure_width,
local_width, nonuniform, function, povm_combo, rgen, eps=1e-10):
"""Verify that estimated probabilities from MPPovmList.est_pmf_from()
are reasonable accurate
.. todo:: This test is too long and should be split into several
smaller tests. Also, some of the testing done here is
redundant.
"""
mps = factory.random_mps(nr_sites, local_dim, rank, rgen)
mps.canonicalize()
sample_povm, fun_povm = povm_combo
estimation_impossible = (sample_povm == "all-y" and
fun_povm in {"local-x", "pauli"})
fromself = sample_povm == fun_povm and measure_width == local_width
# s_povm: POVM used to obtain samples
s_povm = _get_povm(sample_povm, nr_sites, local_dim, measure_width)
# f_povm: POVM on which a linear function is defined
if fromself:
f_povm = s_povm
else:
f_povm = _get_povm(fun_povm, nr_sites, local_dim, local_width)
if function == 'rand':
def coeff(x): return rgen.rand(*x)
elif function == 'randn':
def coeff(x): return rgen.randn(*x)
elif function == 'ones':
def coeff(x): return np.ones(x)
elif function == 'signs':
def coeff(x): return rgen.choice([1., -1.], x)
else:
raise ValueError('Unknown function {!r}'.format(function))
# More POVMs in s_povm means more samples. Consider this in the tests.
n_samples_eff = n_samples * len(s_povm.mpps)
# We divide the coefficients by len(f_povm.mpps) to make the
# estimated value have approximately same magnitude, independently
# of len(f_povm.mpps).
coeff = [coeff(mpp.nsoutdims) / len(f_povm.mpps) for mpp in f_povm.mpps]
samples = tuple(s_povm.sample(
rgen, mps, n_samples, method, mode='mps', eps=eps))
exact_prob = tuple(f_povm.pmf_as_array(mps, 'mps', eps))
# Compute exact estimate directly
exact_est1, exact_var1 = f_povm.lfun([c.ravel() for c in coeff],
None, mps, 'mps', eps)
# Compute exact estimate and variance using the other POVM
exact_est2, exact_var2 = f_povm.lfun_from(s_povm, coeff, mps, 'mps', eps=eps)
if estimation_impossible:
assert np.isnan(exact_est2)
assert np.isnan(exact_var2)
else:
# Estimates must agree.
assert abs(exact_est1 - exact_est2) <= eps
if fromself:
# Variances can be different unless f_povm and s_povm are the same.
assert abs(exact_var1 - exact_var2) <= eps
est, var = f_povm.est_lfun_from(s_povm, coeff, samples, eps)
if fromself:
# In this case, est_lfun() and est_lfun_from() must give exactly
# the same result.
est2, var2 = f_povm.est_lfun([c.ravel() for c in coeff],
None, samples, eps)
assert abs(est - est2) <= eps
assert abs(var - var2) <= eps
# We use est_pmf() to test est_pmf_from()
# again. MPPovmList.est_pmf() just aggregates results from
# MPPovm.est_pmf().
pmf1 = f_povm.est_pmf(samples, normalized=True, eps=eps)
pmf2, _ = zip(*f_povm.est_pmf_from(s_povm, samples, eps=eps))
assert all(abs(p1 - p2).max() <= eps for p1, p2 in zip(pmf1, pmf2))
# The final estimator is based on the samples for
# `s_povm`. Therefore, it is correct to use `n_samples_eff` below
# (and not the "effective samples" for the `f_povm` probability
# estimation returned by :func:`f_povm.est_pmf_from()`).
exact_est = sum(np.inner(c.flat, p.flat) for c, p in zip(coeff, exact_prob))
if estimation_impossible:
assert np.isnan(est)
else:
assert abs(exact_est - exact_est2) <= eps
bound = 20 if fun_povm == 'all-y' else 6
assert abs(est - exact_est) <= bound / n_samples_eff**0.5
if function == 'ones':
assert abs(exact_est - 1) <= eps
assert abs(est - exact_est) <= eps
# The code below will only work for small systems. Probably
# nr_sites = 16 will work, but let's stay safe.
assert nr_sites <= 8, "Larger systems will require a lot of memory"
# Use the estimator from `f_povm._estfun_from_estimator()` to
# compute the exact variance of the estimate. We can assume that
# estimator is mostly correct because we have checked that it
# produces accurate estimates (for large numbers of samples)
# above.
#
# FIXME: Drop the exact variance computation here and use
# exact_var2 from above.
#
# Convert from matching functions + coefficients to coefficients
# for each probability.
n_samples2 = [s.shape[0] for s in samples]
_, est_coeff, est_funs = f_povm._lfun_estimator(s_povm, coeff, n_samples2, eps)
est_p_coeff = [np.zeros(mpp.nsoutdims, float) for mpp in s_povm.mpps]
for fun_coeff, funs, p_coeff, mpp in zip(
est_coeff, est_funs, est_p_coeff, s_povm.mpps):
out = np.unravel_index(range(np.prod(mpp.nsoutdims)), mpp.nsoutdims)
out = np.array(out).T.copy()
for c, fun in zip(fun_coeff, funs):
match = fun(out)
p_coeff.flat[match] += c
exact_prob = tuple(s_povm.pmf_as_array(mps, 'mps', eps))
exact_p_cov = (np.diag(p.flat) - np.outer(p.flat, p.flat) for p in exact_prob)
exact_var = sum(np.inner(c.flat, np.dot(cov, c.flat))
for c, cov in zip(est_p_coeff, exact_p_cov))
if estimation_impossible:
assert np.isnan(var)
else:
assert abs(exact_var - exact_var2) <= eps
if fromself:
# `f_povm` and `s_povm` are equal. We must obtain exactly the
# same result without using the matching functions from above:
exact_prob = tuple(f_povm.pmf_as_array(mps, 'mps', eps))
exact_p_cov = (np.diag(p.flat) - np.outer(p.flat, p.flat)
for p in exact_prob)
exact_var2 = sum(np.inner(c.flat, np.dot(cov, c.flat))
for c, cov in zip(coeff, exact_p_cov))
assert abs(exact_var - exact_var2) <= eps
# Convert variance to variance of the estimator (=average)
exact_var /= n_samples
if sample_povm == 'pauli':
bound = 6
elif fun_povm == 'all-y':
bound = 10
else:
bound = 1
assert n_samples * abs(var - exact_var) <= bound / n_samples_eff**0.5
if function == 'ones':
assert abs(exact_var) <= eps
assert abs(var - exact_var) <= eps
