# API

## Prediction

`BayesianTomography.prediction`

— Method`prediction(outcomes, method::LinearInversion)`

Predict the quantum state from the outcomes of a tomography experiment using the `LinearInversion`

method.

`BayesianTomography.prediction`

— Method```
prediction(outcomes, method::BayesianInference{T};
verbose=false,
σ=T(1e-2),
log_prior=x -> zero(T),
x₀=maximally_mixed_state(Int(√size(method.povm, 2)), T),
nsamples=10^4,
nwarm=10^3,
chain=nothing) where {T}
```

Perform a Bayesian inference on the given `outcomes`

using the `BayesianInference`

`method`

.

**Arguments**

`outcomes`

: The outcomes of the experiment.`method::BayesianInference{T}`

: The Bayesian inference method.`verbose=false`

: Print information about the run.`σ=T(1e-2)`

: The initial standard deviation of the proposal distribution.`log_prior=x -> zero(T)`

: The log-prior function.`x₀=maximally_mixed_state(Int(√size(method.povm, 2)), T)`

: The initial state of the chain.`nsamples=10^4`

: The number of samples to take.`nwarm=10^3`

: The number of warm-up samples to take.`chain=nothing`

: If not`nothing`

, store the chain in this matrix.

**Returns**

A tuple with the mean state, its projection in `method.basis`

and the covariance matrix. The mean state is already returned in matrix form.

`BayesianTomography.LinearInversion`

— Type`LinearInversion(povm, basis=gell_mann_matrices(size(first(povm), 1)))`

Construct a linear inversion method for quantum state tomography.

`BayesianTomography.BayesianInference`

— Type```
BayesianInference(povm::AbstractArray{Matrix{T}},
basis=gell_mann_matrices(size(first(povm), 1), complex(T))) where {T}
```

Create a Bayesian inference object from a POVM.

This is passed to the `prediction`

method in order to perform the Bayesian inference.

## Augmentation

`BayesianTomography.compose_povm`

— Function`compose_povm(povms::AbstractArray{Matrix{T}}...; weights=fill(one(T) / length(povms), length(povms))) where {T}`

Compose a POVM (Positive Operator-Valued Measure) from a set of given POVMs.

**Arguments**

`povms`

: Variable number of POVMs. Each POVM is represented as an array of matrices.`weights`

: An optional array of weights associated with each POVM. If not provided, it defaults to a uniform distribution.

**Returns**

- A new POVM that is a composition of the input POVMs, weighted by their respective weights.

**Example**

```
povm1 = [rand(2,2) for _ in 1:3]
povm2 = [rand(2,2) for _ in 1:3]
composed_povm = compose_povm(povm1, povm2)
```

`BayesianTomography.unitary_transform!`

— Function`unitary_transform!(povm, unitary)`

Apply a unitary transformation to each element of a given POVM (Positive Operator-Valued Measure), modifing it in place.

**Arguments**

`povm`

: The POVM to be transformed. It is represented as an array of matrices.`unitary`

: The unitary matrix representing the transformation to be applied.

**Example**

```
bs_povm = [[1.0+im 0; 0 0], [0 0; 0 1]]
half_wave_plate = [1 1; 1 -1] / √2
unitary_transform!(bs_povm, half_wave_plate)
```

`BayesianTomography.unitary_transform`

— Function`unitary_transform(povm, unitary)`

Apply a unitary transformation to each element of a given POVM (Positive Operator-Valued Measure).

**Arguments**

`povm`

: The POVM to be transformed. It is represented as an array of matrices.`unitary`

: The unitary matrix representing the transformation to be applied.

**Returns**

- A new POVM that is the result of applying the unitary transformation to the input POVM.

**Example**

```
bs_povm = [[1.0+im 0; 0 0], [0 0; 0 1]]
half_wave_plate = [1 1; 1 -1] / √2
unitary_transform!(bs_povm, half_wave_plate)
```

`BayesianTomography.augment_povm`

— Function```
augment_povm(povm::AbstractArray{Matrix{T}}, unitaries...;
weights=fill(one(T) / (length(unitaries) + 1), length(unitaries) + 1) where {T}
```

Augment a POVM (Positive Operator-Valued Measure) by applying a set of unitary transformations to it.

**Arguments**

`povm`

: The POVM to be augmented. It is represented as an array of matrices.`unitaries`

: Variable number of unitary matrices representing the transformations to be applied.`weights`

: An optional array of weights associated with each unitary transformation. If not provided, it defaults to a uniform distribution.

**Returns**

- A new POVM that is the result of applying the unitary transformations to the input POVM.

**Example**

```
bs_povm = [[1.0+im 0; 0 0], [0 0; 0 1]]
half_wave_plate = [1 1; 1 -1] / √2
quater_wave_plate = [1 im; im 1] / √2
povm = augment_povm(bs_povm, half_wave_plate, quater_wave_plate, weights=[1 / 2, 1 / 4, 1 / 4])
```

## Generalized Gell-Mann matrices

`BayesianTomography.triangular_indices`

— Function`triangular_indices(d)`

Generate a vector of tuples representing the indices of the lower triangular part of a square matrix of dimension `d`

.

`BayesianTomography.X_matrix`

— Function`X_matrix(j, k, d, ::Type{T}=ComplexF32) where {T<:Union{Real,Complex}}`

Compute the real off diagonal matrix of the generalized Gell-Mann matrices in dimension `d`

.

The type of the matrix elements is `T`

, which defaults to `ComplexF32`

. The only non-zero elements are `X[j, k] = 1`

and `X[k, j] = 1`

. The matrices are normalized to have unit Hilbert-Schmidt norm.

**Examples**

```
julia> X_matrix(1,2,2)
2×2 Matrix{ComplexF32}:
0.0+0.0im 0.707107+0.0im
0.707107+0.0im 0.0+0.0im
```

`BayesianTomography.Y_matrix`

— Function`Y_matrix(j, k, d, ::Type{T}=ComplexF32) where {T<:Complex}`

Compute the imaginary off diagonal matrix of the generalized Gell-Mann matrices in dimension `d`

.

The type of the matrix elements is `T`

, which defaults to `ComplexF32`

. The only non-zero elements are `Y[j, k] = im`

and `Y[k, j] = -im`

. The matrices are normalized to have unit Hilbert-Schmidt norm.

**Examples**

```
julia> Y_matrix(1,2,2)
2×2 Matrix{ComplexF32}:
0.0+0.0im 0.0+0.707107im
0.0-0.707107im 0.0+0.0im
```

`BayesianTomography.Z_matrix`

— Function`Z_matrix(j, d, ::Type{T}=ComplexF32) where {T<:Union{Real,Complex}}`

Compute the `j`

'th diagonal matrix of the generalized Gell-Mann matrices in dimension `d`

.

The type of the matrix elements is `T`

, which defaults to `ComplexF32`

. The matrices are normalized to have unit Hilbert-Schmidt norm. The identity matrix is returned when `j == 0`

.

**Examples**

```
julia> Z_matrix(0, 3)
3×3 Matrix{ComplexF32}:
0.57735+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.57735+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.57735+0.0im
julia> Z_matrix(1, 3)
3×3 Matrix{ComplexF32}:
0.707107+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im -0.707107+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im 0.0+0.0im
julia> Z_matrix(2, 3)
3×3 Matrix{ComplexF32}:
0.408248+0.0im 0.0+0.0im 0.0+0.0im
0.0+0.0im 0.408248+0.0im 0.0+0.0im
0.0+0.0im 0.0+0.0im -0.816497+0.0im
```

`BayesianTomography.gell_mann_matrices`

— Function`gell_mann_matrices(d, ::Type{T}=ComplexF32; include_identity=true) where {T<:Complex}`

Generate a set of Gell-Mann matrices of dimension `d`

.

The Gell-Mann matrices are a set of `d^2 - 1`

linearly independent, traceless, Hermitian matrices that, when augmented with the identity, form a basis for the space of `d × d`

hermitian matrices.

The matrix order is real off-diagonal (`X_matrix`

), imaginary off-diagonal (`Y_matrix`

) and diagonal (`Z_matrix`

). The off-diagonal matrices follow the order given by `triangular_indices`

.

**Arguments**

`d`

: The dimension of the Gell-Mann matrices.`include_identity`

: A boolean flag indicating whether to include the identity matrix in the set. If this is`true`

, the identity is the first element of the basis

**Returns**

- A 3D array of Gell-Mann matrices. The last dimension is the index of the matrix in the basis.

**Examples**

```
julia> gell_mann_matrices(2,include_identity=false)
2×2×3 Array{ComplexF32, 3}:
[:, :, 1] =
0.0+0.0im 0.707107+0.0im
0.707107+0.0im 0.0+0.0im
[:, :, 2] =
0.0+0.0im 0.0-0.707107im
0.0+0.707107im 0.0+0.0im
[:, :, 3] =
0.707107+0.0im 0.0+0.0im
0.0+0.0im -0.707107+0.0im
```

`BayesianTomography.basis_decomposition`

— Function`basis_decomposition(Ω, basis=gell_mann_matrices(d))`

Decompose the array `Ω`

in the provided orthonormal basis.

If no basis is provided, the Gell-Mann matrices of appropriate dimension are used.

If `Ω`

has dimension d, then `basis`

should be an array with dimesnion `d+1`

with the last dimension indexing the basis elements.

## Representations

`BayesianTomography.History`

— Type`History{T<:Integer}`

A type that represents a history of outcomes.

**Fields**

`history::Vector{T}`

: A vector of outcomes.`history[i]`

is the outcome of the i-th measurement.

`BayesianTomography.reduced_representation`

— Method`reduced_representation(outcomes::Array{T,N}) where {T,N}`

Converts a multi-dimensional array of outcomes into a 2D matrix in which the first row contains the indices of non-zero elements and the second row contains the corresponding non-zero values.

`outcomes`

is a multi-dimensional array of outcomes where the `outcomes[n]`

is the number of times the `n`

-th outcome was observed.

The output is a matrix where the first row contains the indices of non-zero elements from the `outcomes`

array and the second row contains the corresponding non-zero values.

This function has an inverse `complete_representation`

.

**Examples**

```
julia> outcomes = [0, 1, 0, 2, 0, 3]
6-element Vector{Int64}:
0
1
0
2
0
3
julia> reduced_representation(outcomes)
2×3 Matrix{Int64}:
2 4 6
1 2 3
```

`BayesianTomography.reduced_representation`

— Method`reduced_representation(history::History)`

Create a reduced representation of the given history.

**Arguments**

`history::History`

: A History object which contains a history of events.

Return a matrix where each column is a pair (event, count). The event is the unique event from the history and count is the number of times the event has occurred.

**Example**

```
julia> h = History([1,1,1,2,1])
History{Int64}([1, 1, 1, 2, 1])
julia> reduced_representation(h)
2×2 Matrix{Int64}:
2 1
1 4
```

`BayesianTomography.complete_representation`

— Method`complete_representation(outcomes::Matrix{T}, size) where {T}`

Create a complete representation of the given outcomes.

`outcomes`

is a matrix where the first row contains the indices of non-zero elements from of the complete representation and the second row contains the corresponding non-zero values.

Returns a vector of size `size`

where the i-th element is the value of the pair whose index is i in `outcomes`

. If there is no such pair, the value is 0.

This function has an inverse `reduced_representation`

.

**Example**

```
julia> outcomes = [1 2; 3 4]
2×2 Matrix{Int64}:
1 2
3 4
julia> complete_representation(outcomes, (2,2))
2×2 Matrix{Int64}:
3 0
4 0
```

`BayesianTomography.complete_representation`

— Method`complete_representation(history::History{T}, size) where {T}`

Create a complete representation of the given history.

**Arguments**

`history::History`

: A History object which contains a history of outcomes.`size`

: The size of the resulting representation.

**Returns**

`result`

: An array of size`size`

where the i-th element is the number of times the i-th event occurred in the history.

**Example**

```
julia> h = History([1,1,1,2,1])
History{Int64}([1, 1, 1, 2, 1])
julia> complete_representation(h,(2,2))
2×2 Matrix{Int64}:
4 0
1 0
```

## Samplers

`BayesianTomography.sample`

— Function```
sample(type, n_samples)
sample(type)
```

Sample `n_samples`

from `type`

.

If `n_samples`

is not provided, a single sample is returned.

Possible values for type are `HaarUnitary`

, `HaarVector`

, `Simplex`

, `ProductMeasure`

, and `GinibreEnsamble`

.

`BayesianTomography.HaarUnitary`

— Type`HaarUnitary(dim::Int)`

A type representing a Haar-random unitary matrix of dimension `dim`

.

`BayesianTomography.HaarVector`

— Type`HaarVector(dim::Int)`

A type representing a Haar-random unitary vector of dimension `dim`

.

`BayesianTomography.Simplex`

— Type`Simplex(dim::Int)`

A type representing a random point on the simplex embeded in a space of dimension `dim`

.

`BayesianTomography.ProductMeasure`

— Type`ProductMeasure(dim::Int)`

A type representing a measure on the density states. It is a product Haar measure on the unitary group and a uniform (Lebesgue) measure on the simplex.

`BayesianTomography.GinibreEnsamble`

— Type`GinibreEnsamble(dim::Int)`

A type representing a Ginibre ensamble of complex matrices of dimension `dim`

.

## Simulation

`BayesianTomography.simulate_outcomes`

— Function```
simulate_outcomes(ψ::AbstractVector, povm, N; atol=1e-3)
simulate_outcomes(ρ::AbstractMatrix, povm, N; atol=1e-3)
simulate_outcomes(probs, N; atol=1e-3)
```

Simulate the `N`

outcomes of a quantum measurement represented by a `povm`

on a quantum state.

The state can be pure or mixed, and it is represented by a vector `ψ`

or a density matrix `ρ`

, respectively. Alternativelly, one can directly provide the probabilities of the outcomes in the `probs`

array.

`atol`

is the absolute tolerance for the probabilities to be considered non-negative and to sum to 1.

`BayesianTomography.simulate_outcomes!`

— Function`simulate_outcomes!(probs, N; atol=1e-3)`

Simulate the `N`

outcomes of a probability specified by the `probs`

array. The results are stored in the `probs`

array.

`atol`

is the absolute tolerance for the probabilities to be considered non-negative and to sum to 1.

## Misc

`BayesianTomography.fidelity`

— Function```
fidelity(ρ::AbstractMatrix, σ::AbstractMatrix)
fidelity(ψ::AbstractVector, φ::AbstractVector)
```

Calculate the fidelity between two quantum states.

The states can be pure or mixed, and they are represented by vectors `ψ`

and `φ`

or density matrices `ρ`

and `σ`

, respectively.

`BayesianTomography.project2density`

— Function`project2density(ρ)`

Project a Hermitian matrix `ρ`

to a density matrix by setting the negative eigenvalues to zero and normalizing the trace to 1.

`BayesianTomography.project2pure`

— Function`project2pure(ρ)`

Project a Hermitian matrix `ρ`

to a pure state by returning the eigenvector corresponding to the largest eigenvalue.

`BayesianTomography.orthogonal_projection`

— Function`orthogonal_projection(ρ, set)`

Calculate the orthogonal projection of `ρ`

onto `set`

.

`set`

is an array with one more dimension than `ρ`

.

`BayesianTomography.real_orthogonal_projection`

— Function`real_orthogonal_projection(ρ, set)`

Calculate the real part of the orthogonal projection of `ρ`

onto `set`

.

`set`

is an array with one more dimension than `ρ`

.

This function is useful when the projection is expected to be real, but numerical errors may introduce small imaginary parts.

`BayesianTomography.linear_combination`

— Function`linear_combination(xs, set)`

Calculate the linear combination of the elements of `set`

with the coefficients `xs`

.

`BayesianTomography.linear_combination!`

— Function`linear_combination!(ρ, xs, set)`

Calculate the linear combination of the elements of `set`

with the coefficients `xs`

and store the result in `ρ`

.

`LinearAlgebra.isposdef!`

— Method`isposdef!(ρ, xs, set)`

Calculate the linear combination of the elements of `set`

with the coefficients `xs`

and check if the result is a positive definite matrix.

`LinearAlgebra.cond`

— Function`cond(M, p::Real=2)`

Condition number of the matrix `M`

, computed using the operator `p`

-norm. Valid values for `p`

are `1`

, `2`

(default), or `Inf`

.

`cond(povm::Union{AbstractArray{T},AbstractMatrix{T}}, p::Real=2) where {T<:AbstractMatrix}`

Calculate the condition number of the linear transformation associated with the `povm`

.

`BayesianTomography.maximally_mixed_state`

— Function`maximally_mixed_state(d, ::Type{T}) where {T}`

Returns the maximally mixed state of dimension `d`

, represented as a vector of projections in the generalized Gell-Mann basis.

The maximally mixed state is defined as `ρ = I / d`

.

Se also `gell_mann_matrices`

.