# Source code for toqito.random.random_povm

"""Generate random POVM."""
import numpy as np

[docs]def random_povm(dim: int, num_inputs: int, num_outputs: int) -> np.ndarray:
"""
Generate random positive operator valued measurements (POVMs) [WIKPOVM]_.

Examples
==========

We can generate a set of dim-by-dim POVMs consisting of a specific dimension along with
a given number of measurement inputs and measurement outputs. As an example, we can
construct a random set of :math:2-by-:math:2 POVMs of dimension with :math:2 inputs
and :math:2 outputs.

>>> from toqito.random import random_povm
>>> import numpy as np
>>>
>>> dim, num_inputs, num_outputs = 2, 2, 2
>>> povms = random_povm(dim, num_inputs, num_outputs)
>>> povms
[[[[ 0.40313832+0.j,  0.59686168+0.j],
[ 0.91134633+0.j,  0.08865367+0.j]],
[[-0.27285707+0.j,  0.27285707+0.j],
[-0.12086852+0.j,  0.12086852+0.j]]],
[[[-0.27285707+0.j,  0.27285707+0.j],
[-0.12086852+0.j,  0.12086852+0.j]],
[[ 0.452533  +0.j,  0.547467  +0.j],
[ 0.34692158+0.j,  0.65307842+0.j]]]]

We can verify that this constitutes a valid set of POVM elements as checking that these
operators all sum to the identity operator.

>>> np.round(povms[:, :, 0, 0] + povms[:, :, 0, 1])
[[1.+0.j, 0.+0.j],
[0.+0.j, 1.+0.j]]

References
==========
.. [WIKPOVM] Wikipedia: POVM
https://en.wikipedia.org/wiki/POVM

:param dim: The dimensions of the measurements.
:param num_inputs: The number of inputs for the measurement.
:param num_outputs: The number of outputs for the measurement.
:return: A set of dim-by-dim POVMs of shape (dim, dim, num_inputs, num_outputs).
"""
povms = []
gram_vectors = np.random.normal(size=(num_inputs, num_outputs, dim, dim))
for input_block in gram_vectors:
normalizer = sum(
[np.array(output_block).T.conj() @ output_block for output_block in input_block]
)

u_mat, d_mat, _ = np.linalg.svd(normalizer)

output_povms = []
for output_block in input_block:
partial = (
np.array(output_block, dtype=complex).dot(u_mat).dot(np.diag(d_mat ** (-1 / 2.0)))
)
internal = partial.dot(np.diag(np.ones(dim)) ** (1 / 2.0))
output_povms.append(internal.T.conj() @ internal)
povms.append(output_povms)

# This allows us to index the POVMs as [dim, dim, num_inputs, num_outputs].
povms = np.swapaxes(np.array(povms), 0, 2)
povms = np.swapaxes(povms, 1, 3)

return povms