Maths

Nonconvex Phase Synchronization

This is a Python3 implementation of the Nonconvex Phase Synchronisation method found in Boumal [Bou16], page 8.

Methodology and Use Case:

We consider that from measurements, we can only obtain noisy relative phase advances \(\mu_{i} - \mu_{j}\) and want a converging solution to reconstruct the different individual \(\mu_{1}, ...., \mu_{n}\) values.

From measurements, we construct a hermitian matrix C in the shape of:

\[C_{ij} = z_{i} \bar{z}_{j} = e^{i (\mu_{i} - \mu_{j})}\]

A mock one with random values (500 by 500 as we have 500 BPMs per plane in the LHC) would be:

C_matrix = np.exp(1j * np.random.rand(500, 500))

Considering 4 BPMs, the measurement matrix would be:

\[\begin{split}M = \begin{pmatrix} \mu_{1 -> 1} & \mu_{1 -> 2} & \mu_{1 -> 3} & \mu_{1 -> 4} \\ \mu_{2 -> 1} & \mu_{2 -> 2} & \mu_{2 -> 3} & \mu_{2 -> 4} \\ \mu_{3 -> 1} & \mu_{3 -> 2} & \mu_{3 -> 3} & \mu_{3 -> 4} \\ \mu_{4 -> 1} & \mu_{4 -> 2} & \mu_{4 -> 3} & \mu_{4 -> 4} \end{pmatrix}\end{split}\]

Note two particular properties here:

Because our measurements are phase differences, the \(\mathrm{M_matrix}\) will necessarily have zeros on its diagonal (\(\mu_{k -> k} = 0\)).
By definition, since \(\mu_{a -> b} = - \mu_{b -> a}\), \(\mathrm{M_matrix}\) is symmetric.
Also note that for all computations, \(\mathrm{M_matrix}\) needs to be initialised in radians!

We can very simply get our C_matrix (see page 1 of referenced paper) with numpy.exp which, applied to a numpy.ndarray applies the exponential function element-wise.

Then follows:

C_matrix = np.exp(1j * M_matrix)

Note

Since \(\mathrm{M_matrix}\) is symmetric, then C_matrix will be Hermitian. Since \(\mathrm{M_matrix}\) has zeros on its diagonal, C_matrix will have (1 + 0j) on its diagonal.

With added noise to those values (noise should be included in \(\mathrm{M_matrix}\) in the case of measurements), we can reconstruct a good estimator of the original values through the EVM method, provided in the class below.

class pyhdtoolkit.maths.nonconvex_phase_sync.PhaseReconstructor(measurements_hermitian_matrix: numpy.ndarray)[source]

New in version 0.4.0.

Class algorithm to reconstruct your phase values. Make sure to provide vectors as numpy.ndarray with shape (1, N), N being the dimension.

Example

real_noised_measurements = np.ndarray(...)
C_herm = np.exp(1j * np.deg2rad(RealNoisedMeas))
pr = PhaseReconstructor(C_herm)
complex_eigenvector_method_result = pr.get_eigenvector_estimator(
    pr.leading_eigenvector
)
reconstructed = np.abs(
    pr.convert_complex_result_to_phase_values(
        complex_eigenvector_method_result, deg=True
    )
).reshape(n_observation_points)

__init__(measurements_hermitian_matrix: numpy.ndarray) → None[source]

Initialize your reconstructor object from measurements.

Parameters: measurements_hermitian_matrix (np.ndarray) -- a 2D array built from measurements, see module docstring on how to build this matrix.

property alpha: numpy.float64

Factor used to define the new reconstruction matrix. It is taken either as the operator norm of the hermitian noise matrix C_matrix, or as the max value between 0 and the opposite of the min eigenvalue of C_matrix (chosen for implementation, since our noise is included in the measurements). See page 8 of the paper for reference.

Returns: A real scalar value, since C_matrix is Hermitian and the eigenvalues of real symmetric or complex Hermitian matrices are always real (see Strang [Str80], page 222).

c_matrix: numpy.ndarray: Hermitian square matrix from your measurements

c_matrix_eigenvalues: numpy.ndarray: Eigenvalues of c_matrix

c_matrix_eigenvectors: numpy.ndarray: Eigenvectors of c_matrix

static convert_complex_result_to_phase_values(complex_estimator: numpy.ndarray, deg: bool = False) → numpy.ndarray[source]

Casts back the complex form of your result to real phase values.

Parameters

complex_estimator (np.ndarray) -- your reconstructed result’s complex form as a numpy.ndarray.
deg (bool) -- if this is set to True, the result is cast to degrees (from radians) before being returned. Defaults to False.

Returns

A numpy.ndarray with the real phase values of the result.

get_eigenvector_estimator(eigenvector: numpy.ndarray) → numpy.ndarray[source]

Return the eigenvector estimator of a given eigenvector (id est the component-wise projection of said eigenvector onto \(\mathbb{C}_1\). See page 7 of Boumal [Bou16] for the implementation.

Parameters: eigenvector (np.ndarray) -- a numpy array representing the vector.
Returns: A numpy.ndarray object of the same dimension as eigenvector.

property leading_eigenvector: numpy.ndarray

Returns the leading eigenvector of self.c_matrix, which is the eigenvector corresponding to the max eigenvalue (in absolute value).

Returns: A numpy.ndarray object, corresponding to said eigenvector.

reconstruct_complex_phases_evm() → numpy.ndarray[source]

Reconstructs the simplest estimator fom the eigenvector method. The result is in complex form, and will be radians once cast back to real form.

Returns: The complex form of the result as a numpy.ndarray.

property reconstructor_matrix: numpy.ndarray

This is the reconstructor matrix built from self.c_matrix and the alpha property. It is the matrix denoted as \(\widetilde{C}\) on page 8 of Boumal [Bou16].

Returns: A numpy.ndarray, with same dimension as self.c_matrix.

space_dimension: int: Dimension of your measurement space

Stats Fitting

Module implementing methods to find the best fit of statistical distributions to data.

pyhdtoolkit.maths.stats_fitting.best_fit_distribution(data: Union[pandas.core.series.Series, numpy.ndarray], bins: int = 200, ax: matplotlib.axes._axes.Axes = None) → Tuple[scipy.stats._distn_infrastructure.rv_continuous, Tuple[float, ...]][source]

New in version 0.5.0.

Model data by finding the best fit candidate distribution among those in DISTRIBUTIONS. One can find an example use of this function in the gallery.

Parameters

data (Union[pd.Series, np.ndarray]) -- a pandas.Series or numpy.ndarray with your distribution data.
bins (int) -- the number of bins to decompose your data in before fitting.
ax (matplotlib.axes.Axes) -- the matplotlib.axes.Axes on which to plot the probability density function of the different fitted distributions. This should be provided as the axis on which the distribution data is plotted, as it will add to that plot.

Returns

A tuple containing the scipy.stats generator corresponding to the best fit candidate, and the parameters for said generator to best fit the data.

Example

best_fit_func, best_fit_params = best_fit_distribution(data, 200, axis)

pyhdtoolkit.maths.stats_fitting.make_pdf(distribution: scipy.stats._distn_infrastructure.rv_continuous, params: Tuple[float, ...], size: int = 25000) → pandas.core.series.Series[source]

New in version 0.5.0.

Generates a pandas.Series for the distributions’s Probability Distribution Function. This Series will have axis values as index, and PDF values as column. One can find an example use of this function in the gallery.

Parameters

distribution (st.rv_continuous) -- a scipy.stats generator.
params (Tuple[float, ...]) -- the parameters for this generator given back by the fit.
size (int) -- the number of points to evaluate.

Returns

A pandas.Series with the PDF as values, and the corresponding axis values as index.

Example

best_fit_func, best_fit_params = best_fit_distribution(data, 200, axis)
pdf = fitting.make_pdf(best_fit_func, best_fit_params)

pyhdtoolkit.maths.stats_fitting.set_distributions_dict(dist_dict: Dict[scipy.stats._distn_infrastructure.rv_continuous, str]) → None[source]

New in version 0.5.0.

Sets DISTRIBUTIONS as the provided dict. This allows the user to define the distributions to try and fit against the data. One can find an example use of this function in the gallery.

Parameters: dist_dict (Dict[st.rv_continuous, str]) -- dictionnary with the wanted distributions, in the format of DISTRIBUTIONS, aka with scipy.stats generator objects as keys, and a string representation of their name as value.
Returns: Nothing, but modifies the global DISTRIBUTIONS dict called by other functions.

Example

import scipy.stats as st
tested_dists = {st.chi: "Chi", st.expon: "Exponential", st.laplace: "Laplace"}
set_distributions_dict(tested_dists)

Utilities

Module with utility functions used throughout the nonconvex_phase_sync and stats_fitting modules.

pyhdtoolkit.maths.utils.get_magnitude(value: float) → int[source]

New in version 0.8.2.

Returns the determined magnitude of the provided value. This corresponds to the power of 10 that would be necessary to reduce value to a \(X \cdot 10^{n}\) form. In this case, n is the result.

Parameters: value (float) -- value to determine the magnitude of.
Returns: The magnitude of the provided value, as an int.

Examples

get_magnitude(10)
# returns 1

get_magnitude(0.0311)
# returns -2

get_magnitude(1e-7)
# returns -7

pyhdtoolkit.maths.utils.get_scaled_values_and_magnitude_string(values_array: Union[pandas.core.frame.DataFrame, numpy.ndarray], force_magnitude: float = None) → Tuple[Union[pandas.core.frame.DataFrame, numpy.ndarray], str][source]

New in version 0.8.2.

Conveniently scales the provided values to the best determined magnitude, and returns the scaled values and the magnitude string to use in plots labels.

Parameters

values_array (Union[pd.DataFrame, np.ndarray]) -- vectorised structure containing the values to scale.
force_magnitude (float) -- a specific magnitude value to use for the scaling, if desired.

Returns

A tuple of the scaled values (same type as the provided ones) and the string to use for the scale in plots labels and legends.

Example

import numpy as np
q = np.array([-330,  230,  430, -720,  750, -110,  410, -340, -950, -630])
get_scaled_values_and_magnitude_string(q)
# returns (array([-3.3,  2.3,  4.3, -7.2,  7.5, -1.1,  4.1, -3.4, -9.5, -6.3]), '{-2}')