Source code for PyDynamic.misc.tools

"""A collection of miscellaneous helper functions

This module contains the following functions:

* :func:`FreqResp2RealImag`: Calculate real and imaginary parts from frequency
  response
* :func:`is_2d_matrix`: Check if a np.ndarray is a matrix
* :func:`is_2d_square_matrix`: Check if a np.ndarray is a two-dimensional square matrix
* :func:`is_vector`: Check if a np.ndarray is a vector
* :func:`make_semiposdef`: Make quadratic matrix positive semi-definite
* :func:`normalize_vector_or_matrix`: Scale an array of numbers to the interval between
  zero and one
* :func:`number_of_rows_equals_vector_dim`: Check if a matrix and a vector match in size
* :func:`plot_vectors_and_covariances_comparison`: Plot two vectors and their
  covariances side-by-side for visual comparison
* :func:`print_mat`: Print matrix (2D array) to the console or return as formatted
  string
* :func:`print_vec`: Print vector (1D array) to the console or return as formatted
  string
* :func:`progress_bar`: A simple and reusable progress-bar
* :func:`shift_uncertainty`: Shift the elements in the vector x and associated
  uncertainties ux
* :func:`trimOrPad`: trim or pad (with zeros) a vector to desired length
* :func:`complex_2_real_imag`: Take a np.ndarray with dtype complex and return
  real and imaginary parts
* :func:`real_imag_2_complex`: Take a np.ndarray with real and imaginary parts
  and return dtype complex ndarray
* :func:`separate_real_imag_of_mc_samples`: Split a np.ndarray containing MonteCarlo
  samples' real and imaginary parts
* :func:`separate_real_imag_of_vector`: Split a np.ndarray containing real and
  imaginary parts into half
"""

__all__ = [
    "print_mat",
    "print_vec",
    "make_semiposdef",
    "FreqResp2RealImag",
    "make_equidistant",
    "trimOrPad",
    "progress_bar",
    "shift_uncertainty",
    "is_vector",
    "is_2d_matrix",
    "number_of_rows_equals_vector_dim",
    "plot_vectors_and_covariances_comparison",
    "is_2d_square_matrix",
    "normalize_vector_or_matrix",
    "complex_2_real_imag",
    "real_imag_2_complex",
    "separate_real_imag_of_mc_samples",
    "separate_real_imag_of_vector",
]

from typing import Any, List, Optional, Union

import numpy as np
from matplotlib import pyplot as plt
from matplotlib.colors import Normalize
from scipy.sparse import eye, issparse
from scipy.sparse.linalg import eigs


[docs]def shift_uncertainty(x: np.ndarray, ux: np.ndarray, shift: int): """Shift the elements in the vector x and associated uncertainties ux This function uses :func:`numpy.roll` to shift the elements in x and ux. See the linked official documentation for details. Parameters ---------- x : np.ndarray of shape (N,) vector of estimates ux : float, np.ndarray of shape (N,) or of shape (N,N) uncertainty associated with the vector of estimates shift : int amount of shift Returns ------- shifted_x : (N,) np.ndarray shifted vector of estimates shifted_ux : float, np.ndarray of shape (N,) or of shape (N,N) uncertainty associated with the shifted vector of estimates Raises ------ ValueError If shift, x or ux are of unexpected type, dimensions of x and ux do not fit or ux is of unexpected shape """ shifted_x = _shift_vector(vector=x, shift=shift) if isinstance(ux, float): return shifted_x, ux if isinstance(ux, np.ndarray): if is_vector(ux): return shifted_x, _shift_vector(vector=ux, shift=shift) elif is_2d_square_matrix(ux): shifted_ux = _shift_2d_matrix(ux, shift) return shifted_x, shifted_ux raise ValueError( "shift_uncertainty: input uncertainty ux is expected to be a vector or " f"two-dimensional, square matrix but is of shape {ux.shape}." ) raise ValueError( "shift_uncertainty: input uncertainty ux is expected to be a float or a " f"numpy.ndarray but is of type {type(ux)}." )
def _cast_shift_to_int(shift: Any) -> int: try: return int(shift) except ValueError: raise ValueError( "shift_uncertainty: shift is expected to be type int or at least " f"cast-able to int, but is {shift} of type {type(shift)}. Please provide " "a valid value." ) def _shift_vector(vector: np.ndarray, shift: int) -> np.ndarray: return np.roll(vector, shift) def _shift_2d_matrix(matrix: np.ndarray, shift: int) -> np.ndarray: return np.roll(matrix, (shift, shift), axis=(0, 1))
[docs]def trimOrPad( array: Union[List, np.ndarray], length: int, mode: Optional[str] = "constant" ): """Trim or pad (with zeros) a vector to the desired length Either trim or zero-pad an array to achieve the required `length`. Both actions are applied to the end of the array. Parameters ---------- array : list, 1D np.ndarray original data length : int length of output mode : str, optional handed over to np.pad, default "constant" Returns ------- array_modified : np.ndarray of shape (length,) An either trimmed or zero-padded array """ if len(array) < length: # pad zeros to the right if too short return np.pad(array, (0, length - len(array)), mode=mode) else: # trim to given length otherwise return array[0:length]
[docs]def make_semiposdef( matrix: np.ndarray, maxiter: Optional[int] = 10, tol: Optional[float] = 1e-12, verbose: Optional[bool] = False, ) -> np.ndarray: """Make quadratic matrix positive semi-definite by increasing its eigenvalues Parameters ---------- matrix : array_like of shape (N,N) the matrix to process maxiter : int, optional the maximum number of iterations for increasing the eigenvalues, defaults to 10 tol : float, optional tolerance for deciding if pos. semi-def., defaults to 1e-12 verbose : bool, optional If True print smallest eigenvalue of the resulting matrix, if False (default) be quiet Returns ------- (N,N) array_like quadratic positive semi-definite matrix Raises ------ ValueError If matrix is not square. """ n, m = matrix.shape if n != m: raise ValueError("Matrix has to be quadratic") # use specialised functions for sparse matrices if issparse(matrix): # enforce symmetric matrix matrix = 0.5 * (matrix + matrix.T) # calculate smallest eigenvalue e = np.min(np.real(eigs(matrix, which="SR", return_eigenvectors=False))) count = 0 # increase the eigenvalues until matrix is positive semi-definite while e < tol and count < maxiter: matrix += (np.absolute(e) + tol) * eye(n, format=matrix.format) e = np.min(np.real(eigs(matrix, which="SR", return_eigenvectors=False))) count += 1 e = np.min(np.real(eigs(matrix, which="SR", return_eigenvectors=False))) # same procedure for non-sparse matrices else: matrix = 0.5 * (matrix + matrix.T) count = 0 e = np.min(np.real(np.linalg.eigvals(matrix))) while e < tol and count < maxiter: e = np.min(np.real(np.linalg.eigvals(matrix))) matrix += (np.absolute(e) + tol) * np.eye(n) e = np.min(np.real(np.linalg.eigvals(matrix))) if verbose: print("Final result of make_semiposdef: smallest eigenvalue is %e" % e) return matrix
[docs]def FreqResp2RealImag( Abs: np.ndarray, Phase: np.ndarray, Unc: np.ndarray, MCruns: Optional[int] = 1000 ): """Calculate real and imaginary parts from frequency response Calculate real and imaginary parts from amplitude and phase with associated uncertainties. Parameters ---------- Abs : (N,) array_like amplitude values Phase : (N,) array_like phase values in rad Unc : (2N, 2N) or (2N,) array_like uncertainties either as full covariance matrix or as its main diagonal MCruns : int, optional number of iterations for Monte Carlo simulation, defaults to 1000 Returns ------- Re, Im : (N,) array_like best estimate of real and imaginary parts URI : (2N, 2N) array_like uncertainties assoc. with Re and Im """ if len(Abs) != len(Phase) or 2 * len(Abs) != len(Unc): raise ValueError("\nLength of inputs are inconsistent.") if len(Unc.shape) == 1: Unc = np.diag(Unc) Nf = len(Abs) AbsPhas = np.random.multivariate_normal( np.hstack((Abs, Phase)), Unc, int(MCruns) ) # draw MC inputs H = AbsPhas[:, :Nf] * np.exp( 1j * AbsPhas[:, Nf:] ) # calculate complex frequency response values RI = np.hstack((np.real(H), np.imag(H))) # transform to real, imag Re = np.mean(RI[:, :Nf]) Im = np.mean(RI[:, Nf:]) URI = np.cov(RI, rowvar=False) return Re, Im, URI
[docs]def make_equidistant(*args, **kwargs): """ .. deprecated:: 2.0.0 Please use :func:`PyDynamic.uncertainty.interpolate.make_equidistant` """ raise DeprecationWarning( "The method `PyDynamic.misc.tools.make_equidistant` is moved " "to :mod:`PyDynamic.uncertainty.interpolate.make_equidistant` since the last " "major release 2.0.0. Please switch to the current module immediately and use " "the current function " ":func:`PyDynamic.uncertainty.interpolate.make_equidistant`. Please change " "'from PyDynamic.misc.tools import make_equidistant' to 'from " "PyDynamic.uncertainty.interpolate import make_equidistant'." )
[docs]def progress_bar( count, count_max, width: Optional[int] = 30, prefix: Optional[str] = "", done_indicator: Optional[str] = "#", todo_indicator: Optional[str] = ".", fout: Optional = None, ): """A simple and reusable progress-bar Parameters ---------- count : int current status of iterations, assumed to be zero-based count_max : int total number of iterations width : int, optional width of the actual progressbar (actual printed line will be wider), default to 30 prefix : str, optional some text that will be printed in front of the bar (i.e. "Progress of ABC:"), if not given only progressbar itself will be printed done_indicator : str, optional what character is used as "already-done"-indicator, defaults to "#" todo_indicator : str, optional what character is used as "not-done-yet"-indicator, defaults to "." fout : file-object, optional where the progress-bar should be written/printed to, defaults to direct print to stdout """ x = int(width * (count + 1) / count_max) progressString = "{PREFIX}[{DONE}{NOTDONE}] {COUNT}/{COUNTMAX}\r".format( PREFIX=prefix, DONE=x * done_indicator, NOTDONE=(width - x) * todo_indicator, COUNT=count + 1, COUNTMAX=count_max, ) if fout is not None: fout.write(progressString) else: print(progressString)
[docs]def is_vector(ndarray: np.ndarray) -> bool: """Check if a np.ndarray is a vector, i.e. is of shape (n,) Parameters ---------- ndarray : np.ndarray the array to check Returns ------- bool True, if the array expands over one dimension only, False otherwise """ return len(ndarray.shape) == 1
[docs]def is_2d_matrix(ndarray: np.ndarray) -> bool: """Check if a np.ndarray is a matrix, i.e. is of shape (n,m) Parameters ---------- ndarray : np.ndarray the array to check Returns ------- bool True, if the array expands over exactly two dimensions, False otherwise """ return len(ndarray.shape) == 2
[docs]def number_of_rows_equals_vector_dim(matrix: np.ndarray, vector: np.ndarray) -> bool: """Check if a matrix has the same number of rows as a vector Parameters ---------- matrix : np.ndarray the matrix, that is supposed to have the same number of rows vector : np.ndarray the vector, that is supposed to have the same number of elements Returns ------- bool True, if the number of rows coincide, False otherwise """ return len(vector) == matrix.shape[0]
[docs]def plot_vectors_and_covariances_comparison( vector_1: np.ndarray, vector_2: np.ndarray, covariance_1: np.ndarray, covariance_2: np.ndarray, title: Optional[str] = "Comparison between two vectors and corresponding " "uncertainties", label_1: Optional[str] = "vector_1", label_2: Optional[str] = "vector_2", ): """Plot two vectors and their covariances side-by-side for visual comparison Parameters ---------- vector_1 : np.ndarray the first vector to compare vector_2 : np.ndarray the second vector to compare covariance_1 : np.ndarray the first covariance matrix to compare covariance_2 : np.ndarray the second covariance matrix to compare title : str, optional the title for the comparison plot, defaults to `"Comparison between two vectors and corresponding uncertainties"` label_1 : str, optional the label for the first vector in the legend and title for the first covariance plot, defaults to "vector_1" label_2 : str, optional the label for the second vector in the legend and title for the second covariance plot, defaults to "vector_2" """ fig, ax = plt.subplots(nrows=2, ncols=2) fig.suptitle(title) ax[0][0].imshow(covariance_1) ax[0][0].set_title(label_1 + " uncertainties") ax[0][1].imshow(covariance_2) ax[0][1].set_title(label_2 + " uncertainties") ax[1][0].plot(vector_1, label=label_1) ax[1][0].plot(vector_2, label=label_2) ax[1][0].legend() ax[1][0].set_title(label_1 + " and " + label_2) ax[1][1].imshow(covariance_2 - covariance_1, norm=Normalize()) ax[1][1].set_title("Relative difference of uncertainties") plt.show()
[docs]def is_2d_square_matrix(ndarray: np.ndarray) -> bool: """Check if a np.ndarray is a two-dimensional square matrix, i.e. is of shape (n,n) Parameters ---------- ndarray : np.ndarray the array to check Returns ------- bool True, if the array expands over exactly two dimensions of similar size, False otherwise """ return is_2d_matrix(ndarray) and ndarray.shape[0] == ndarray.shape[1]
[docs]def normalize_vector_or_matrix(numbers: np.ndarray) -> np.ndarray: """Scale an array of numbers to the interval between zero and one If all values in the array are the same, the output array will be constant zero. Parameters ---------- numbers : np.ndarray the :class:`numpy.ndarray` to normalize Returns ------- np.ndarray the normalized array """ minimum = translator = np.min(numbers) array_span = np.max(numbers) - minimum normalizer = array_span or 1.0 return (numbers - translator) / normalizer
[docs]def complex_2_real_imag(array: np.ndarray) -> np.ndarray: r"""Take an array of any non-flexible scalar dtype to return real and imaginary part The input array :math:`x \in \mathbb R^m` is reassembled to the form of the expected input of some of the functions in the modules :mod:`propagate_DFT <PyDynamic.uncertainty.propagate_DFT>` and :mod:`fit_filter <PyDynamic.model_estimation.fit_filter>`: :math:`y = \left( \operatorname{Re}(x), \operatorname{Im}(x) \right)`. Parameters ---------- array : np.ndarray of shape (M,) the array to assemble the version with real and imaginary parts from Returns ------- np.ndarray of shape (2M,) the array of real and imaginary parts """ return np.hstack((np.real(array), np.imag(array)))
[docs]def real_imag_2_complex(array: np.ndarray) -> np.ndarray: r"""Take a np.ndarray with real and imaginary parts and return dtype complex ndarray The input array :math:`x \in \mathbb R^{2m}` representing a complex vector :math:`y \in \mathbb C^m` has the form of the expected input of some of the functions in the modules :mod:`propagate_DFT <PyDynamic.uncertainty.propagate_DFT>` and :mod:`fit_filter <PyDynamic.model_estimation.fit_filter>`: :math:`x = \left( \operatorname{Re}(y), \operatorname{Im}(y) \right)` or a np.ndarray containing several of these. Parameters ---------- array : np.ndarray of shape (N,2M) or of shape (2M,) the array of any integer or floating dtype to assemble the complex version of Returns ------- np.ndarray of shape (N,M) or of shape (M,) the complex array """ if is_2d_matrix(array): real, imag = separate_real_imag_of_mc_samples(array) else: real, imag = separate_real_imag_of_vector(array) return real + 1j * imag
[docs]def separate_real_imag_of_mc_samples(array: np.ndarray) -> List[np.ndarray]: r"""Split a np.ndarray containing MonteCarlo samples real and imaginary parts The input array :math:`x \in \mathbb R^{n \times 2m}` representing an n-elemental array of complex vectors :math:`y_i \in \mathbb C^m` has the form of the expected input of some of the functions in the modules :mod:`propagate_DFT <PyDynamic.uncertainty.propagate_DFT>` and :mod:`fit_filter <PyDynamic.model_estimation.fit_filter>`: :math:`x = \left( \operatorname{Re}(y_i), \operatorname{Im}(y_i) \right)_{i=1,\ldots,n}`. Parameters ---------- array : np.ndarray of shape (N,2M) the array of any integer or floating dtype to assemble the complex version of Returns ------- list of two np.ndarrays of shape (N,M) two-element list of the two arrays containing the real and imaginary parts """ if _vector_has_odd_length(array[0]): raise ValueError( "separate_real_imag_of_mc_samples: vectors of real and imaginary " "parts are expected to contain exactly as many real as " f"imaginary parts but the first one is of odd length={len(array[0])}." ) return np.split(ary=array, indices_or_sections=2, axis=1)
[docs]def separate_real_imag_of_vector(vector: np.ndarray) -> List[np.ndarray]: r"""Split a np.ndarray containing real and imaginary parts into half The input array :math:`x \in \mathbb R^{2m}` representing a complex vector :math:`y \in \mathbb C^m` has the form of the expected input of some of the functions in the modules :mod:`propagate_DFT <PyDynamic.uncertainty.propagate_DFT>` and :mod:`fit_filter <PyDynamic.model_estimation.fit_filter>`: :math:`x = \left( \operatorname{Re}(y), \operatorname{Im}(y) \right)`. Parameters ---------- vector : np.ndarray of shape (2M,) the array of any integer or floating dtype to assemble the complex version of Returns ------- list of two np.ndarrays of shape (M,) two-element list of the two arrays containing the real and imaginary parts """ if _vector_has_odd_length(vector): raise ValueError( "separate_real_imag_of_vector: vector of real and imaginary " "parts is expected to contain exactly as many real as " f"imaginary parts but is of odd length={len(vector)}." ) return np.split(ary=vector, indices_or_sections=2)
def _vector_has_odd_length(vector: np.ndarray) -> bool: return len(vector) % 2 == 1