15472938725

from copy import deepcopy
from typing import Any

import numpy as np
import numpy.typing as npt
import scipy


def get_rotation_matrix(
    vec_start: npt.NDArray[np.float64], vec_end: npt.NDArray[np.float64]
) -> npt.NDArray[np.float64]:
    """
    Calculate the rotation matrix to align two unit vectors.

    Given a two (unit) vectors, vec_start and vec_end, this function calculates the
    rotation matrix U, so that U * vec_start = vec_end.

    U the is rotation matrix that rotates vec_start to point in the direction
    of vec_end.

    https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/897677

    Parameters
    ----------
    vec_start, vec_end : NDArray[float64]
        Two vectors that should be aligned. Both vectors must have a l2-norm of 1.

    Returns
    -------
    NDArray[float64]
        The rotation matrix U as npt.NDArray with shape (3,3)
    """
    if not np.isclose(np.linalg.norm(vec_start), 1) and not np.isclose(
        np.linalg.norm(vec_end), 1
    ):
        raise ValueError("`vec_start` and `vec_end` args must be unit vectors")

    v = np.cross(vec_start, vec_end)
    c = np.dot(vec_start, vec_end)
    v_x = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])
    return np.eye(3) + v_x + v_x.dot(v_x) / (1 + c)


def get_rotation_matrix_around_axis(axis: npt.NDArray, phi: float) -> npt.NDArray:
    """
    Generate a rotation matrix around a given vector.

    Parameters
    ----------
    axis : NDArray
        Axis around which the rotation is done.
    phi : float
        Angle of rotation around axis in radiants.

    Returns
    -------
    NDArray
        Rotation matrix
    """
    axis_vec = np.array(axis, dtype=np.float64)
    axis_vec /= np.linalg.norm(axis_vec)

    eye = np.eye(3, dtype=np.float64)
    ddt = np.outer(axis_vec, axis_vec)
    skew = np.array(
        [
            [0, axis_vec[2], -axis_vec[1]],
            [-axis_vec[2], 0, axis_vec[0]],
            [axis_vec[1], -axis_vec[0], 0],
        ],
        dtype=np.float64,
    )

    return ddt + np.cos(phi) * (eye - ddt) + np.sin(phi) * skew


def get_rotation_matrix_around_z_axis(phi: float) -> npt.NDArray:
    """
    Generate a rotation matrix around the z axis.

    Parameters
    ----------
    phi : float
        Angle of rotation around axis in radiants.

    Returns
    -------
    NDArray
        Rotation matrix
    """
    return get_rotation_matrix_around_axis(np.array([0.0, 0.0, 1.0]), phi)


def get_mirror_matrix(normal_vector: npt.NDArray) -> npt.NDArray:
    """
    Generate a transformation matrix for mirroring through plane given by the normal.

    Parameters
    ----------
    normal_vector : NDArray
        Normal vector of the mirror plane.

    Returns
    -------
    NDArray
        Mirror matrix
    """
    n_vec = normal_vector / np.linalg.norm(normal_vector)
    eps = np.finfo(np.float64).eps
    a = n_vec[0]
    b = n_vec[1]
    c = n_vec[2]
    M = np.array(
        [
            [1 - 2 * a**2, -2 * a * b, -2 * a * c],
            [-2 * a * b, 1 - 2 * b**2, -2 * b * c],
            [-2 * a * c, -2 * b * c, 1 - 2 * c**2],
        ]
    )
    M[np.abs(M) < eps * 10] = 0
    return M


def get_angle_between_vectors(
    vector_1: npt.NDArray, vector_2: npt.NDArray
) -> npt.NDArray:
    """
    Determine angle between two vectors.

    Parameters
    ----------
    vector_1 : NDArray
    vector_2 : NDArray

    Returns
    -------
    float
        Angle in radiants.
    """
    return (
        np.dot(vector_1, vector_2) / np.linalg.norm(vector_1) / np.linalg.norm(vector_2)
    )


def get_fractional_coords(
    cartesian_coords: npt.NDArray, lattice_vectors: npt.NDArray
) -> npt.NDArray:
    """
    Transform cartesian coordinates into fractional coordinates.

    Parameters
    ----------
    cartesian_coords: NDArray
        Cartesian coordinates of atoms (can be Nx2 or Nx3)
    lattice_vectors: [N_dim x N_dim] numpy array:
        Matrix of lattice vectors: Each ROW corresponds to one lattice vector!

    Returns
    -------
    fractional_coords: [N x N_dim] numpy array
        Fractional coordinates of atoms
    """
    fractional_coords = np.linalg.solve(lattice_vectors.T, cartesian_coords.T)
    return fractional_coords.T


def get_cartesian_coords(
    frac_coords: npt.NDArray, lattice_vectors: npt.NDArray
) -> npt.NDArray:
    """
    Transform fractional coordinates into cartesian coordinates.

    Parameters
    ----------
    frac_coords: NDArray
        Fractional coordinates of atoms (can be Nx2 or Nx3)
    lattice_vectors: NDArray
        Matrix of lattice vectors: Each ROW corresponds to one lattice vector!

    Returns
    -------
    NDArray
        Cartesian coordinates of atoms

    """
    return np.dot(frac_coords, lattice_vectors)


def get_triple_product(a: npt.NDArray, b: npt.NDArray, c: npt.NDArray) -> npt.NDArray:
    """
    Get the triple product (DE: Spatprodukt): a*(bxc).

    Parameters
    ----------
    a: NDArray
        TODO
    b: NDArray
        TODO
    c: NDArray
        TODO

    Returns
    -------
    NDarray
        Triple product of each input vector
    """
    if len(a) != 3 or len(b) != 3 or len(c) != 3:
        raise ValueError("Each vector must be of length 3")

    return np.dot(np.cross(a, b), c)


def smooth_function(y: npt.NDArray, box_pts: int) -> npt.NDArray[np.floating[Any]]:
    """
    Smooths a function using convolution.

    Parameters
    ----------
    y : NDArray
        TODO
    box_pts : int
        TODO

    Returns
    -------
    y_smooth : NDArray[floating[Any]]
        TODO
    """
    box = np.ones(box_pts) / box_pts
    return np.convolve(y, box, mode="same")


def get_cross_correlation_function(
    signal_0: npt.NDArray,
    signal_1: npt.NDArray,
    detrend: bool = False,
) -> npt.NDArray:
    """
    Calculate the autocorrelation function for a given signal.

    Parameters
    ----------
    signal_0 : NDArray
        First siganl for which the correlation function should be calculated.
    signal_1 : NDArray
        Second siganl for which the correlation function should be calculated.

    Returns
    -------
    NDArray
        Autocorrelation function from 0 to max_lag.
    """
    if detrend:
        signal_0 = scipy.signal.detrend(signal_0)
        signal_1 = scipy.signal.detrend(signal_1)

    cross_correlation = np.correlate(signal_0, signal_1, mode="full")
    cross_correlation = cross_correlation[cross_correlation.size // 2 :]

    # normalize by number of overlapping data points
    cross_correlation /= np.arange(cross_correlation.size, 0, -1)
    cutoff = int(cross_correlation.size * 0.75)
    return cross_correlation[:cutoff]


def get_autocorrelation_function_manual_lag(
    signal: npt.NDArray, max_lag: int
) -> npt.NDArray:
    """
    Alternative method for autocorrelation of a signal using numpy.corrcoef.

    Allows the lag to be set manually.

    Parameters
    ----------
    signal : NDArray
        Siganl for which the autocorrelation function should be calculated.
    max_lag : int | None
        Autocorrelation will be calculated for a range of 0 to max_lag,
        where max_lag is the largest lag for the calculation of the
        autocorrelation function

    Returns
    -------
    autocorrelation : npt.NDArray
        Autocorrelation function from 0 to max_lag.
    """
    lag = npt.NDArray(range(max_lag))
    autocorrelation = np.array([np.nan] * max_lag)

    for i in lag:
        corr = 1.0 if i == 0 else np.corrcoef(signal[i], signal[:-i])[0][1]
        autocorrelation[i] = corr

    return autocorrelation


def get_fourier_transform(
    signal: npt.NDArray, time_step: float
) -> tuple[npt.NDArray, npt.NDArray]:
    """
    Calculate the fourier transform of a given siganl.

    Parameters
    ----------
    signal : NDArray
        Siganl for which the autocorrelation function should be calculated.
    time_step : float
        Time step of the signal in seconds.

    Returns
    -------
    tuple[NDArray, NDArray]
        Frequencs and absolute values of the fourier transform.

    """
    f = scipy.fft.fftfreq(signal.size, d=time_step)
    y = scipy.fft.fft(signal)

    L = f >= 0

    return f[L], y[L]


def lorentzian(
    x: float | npt.NDArray[np.float64], a: float, b: float, c: float
) -> float | npt.NDArray[np.float64]:
    """
    Return a Lorentzian function.

    Parameters
    ----------
    x : float | NDArray[float64]
        Argument x of f(x) --> y.
    a : float
        Maximum of Lorentzian.
    b : float
        Width of Lorentzian.
    c : float
        Magnitude of Lorentzian.

    Returns
    -------
    f : float | NDArray[float64]
        Outupt of a Lorentzian function.
    """
    return c / (1.0 + ((x - a) / (b / 2.0)) ** 2)


def exponential(
    x: float | npt.NDArray[np.float64], a: float, b: float
) -> float | npt.NDArray[np.float64]:
    """
    Return an exponential function.

    Parameters
    ----------
    x : float | NDArray[float64]
        Argument x of f(x) --> y
    a : float
        TODO
    b : float
        TODO
    """
    return a * np.exp(x * b)


def double_exponential(
    x: float | npt.NDArray, a: float, b: float, c: float
) -> float | npt.NDArray:
    """TODO."""
    return a * (np.exp(x * b) + np.exp(x * c))


def gaussian_window(N: int, std: float = 0.4) -> npt.NDArray[np.float64]:
    """
    Generate a Gaussian window.

    Parameters
    ----------
    N : int
        Number of points in the window.
    std : float
        Standard deviation of the Gaussian window, normalized
        such that the maximum value occurs at the center of the window.

    Returns
    -------
    window : NDarray[float64]
        Gaussian window of length N.

    """
    n = np.linspace(-1, 1, N)
    return np.exp(-0.5 * (n / std) ** 2)


def apply_gaussian_window(
    data: npt.NDArray[np.float64], std: float = 0.4
) -> npt.NDArray[np.float64]:
    """
    Apply a Gaussian window to an array.

    Parameters
    ----------
    data : NDarray[float64]
        Input data array to be windowed.
    std : float
        Standard deviation of the Gaussian window.

    Returns
    -------
    windowed_data : NDArray[float64]
        Windowed data array.
    """
    N = len(data)
    window = gaussian_window(N, std)
    return data * window


def hann_window(N: int) -> npt.NDArray[np.float64]:
    """
    Generate a Hann window.

    Parameters
    ----------
    N : int
        Number of points in the window.

    Returns
    -------
    NDArray
        Hann window of length N.
    """
    return 0.5 * (1 - np.cos(2 * np.pi * np.arange(N) / (N - 1)))


def apply_window(data: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:
    """
    Apply a Hann window to an array.

    Parameters
    ----------
    data : NDarray[float64]
        Input data array to be windowed.

    Returns
    -------
    NDarray[float64]
        Windowed data array.
    """
    N = len(data)
    window = hann_window(N)
    return data * window


def norm_matrix_by_dagonal(matrix: npt.NDArray) -> npt.NDArray:
    """
    Norms a matrix such that the diagonal becomes 1.

    | a_11 a_12 a_13 |       |   1   a'_12 a'_13 |
    | a_21 a_22 a_23 |  -->  | a'_21   1   a'_23 |
    | a_31 a_32 a_33 |       | a'_31 a'_32   1   |

    Parameters
    ----------
    matrix : NDArray
        Matrix that should be normed.

    Returns
    -------
    matrix : NDArray
        Normed matrix.
    """
    diagonal = np.array(np.diagonal(matrix))
    L = diagonal == 0.0
    diagonal[L] = 1.0

    new_matrix = deepcopy(matrix)
    new_matrix /= np.sqrt(
        np.tile(diagonal, (matrix.shape[1], 1)).T
        * np.tile(diagonal, (matrix.shape[0], 1))
    )

    return new_matrix


def mae(delta: npt.NDArray) -> np.floating:
    """
    Calculate the mean absolute error from a list of value differnces.

    Parameters
    ----------
    delta : NDArray
        Array containing differences

    Returns
    -------
    float
        mean absolute error
    """
    return np.mean(np.abs(delta))


def rel_mae(
    delta: npt.NDArray, target_val: npt.NDArray
) -> np.floating | np.float64 | float:
    """
    Compute relative MAE from value differences and corresponding target values.

    Parameters
    ----------
    delta : NDArray
        Array containing differences
    target_val : NDArray
        Array of target values against which the difference should be compared

    Returns
    -------
    float
        relative mean absolute error
    """
    target_norm = np.mean(np.abs(target_val))
    return np.mean(np.abs(delta)).item() / (target_norm + 1e-9)


def rmse(delta: npt.NDArray) -> float:
    """
    Calculate the root mean sqare error from a list of value differnces.

    Parameters
    ----------
    delta : NDArray
        Array containing differences

    Returns
    -------
    float
        root mean square error
    """
    return np.sqrt(np.mean(np.square(delta)))


def rel_rmse(delta: npt.NDArray, target_val: npt.NDArray) -> float:
    """
    Calculate the relative root mean sqare error from a list of value differences.

    Parameters
    ----------
    delta : NDArray
        Array containing differences
    target_val : NDArray
        Array of target values against which the difference should be compared

    Returns
    -------
    float
        relative root mean sqare error

    """
    target_norm = np.sqrt(np.mean(np.square(target_val)))
    return np.sqrt(np.mean(np.square(delta))) / (target_norm + 1e-9)


def get_moving_average(
    signal: npt.NDArray[np.float64], window_size: int
) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
    """
    Cacluated the moving average and the variance around the moving average.

    Parameters
    ----------
    signal : NDArray[float64]
        Signal for which the moving average should be calculated.
    window_size : int
        Window size for the mocing average.

    Returns
    -------
    npt.NDArray[floating]
        Moving average.
    npt.NDArray[floating]
        Variance around the moving average.
    """
    moving_avg = np.convolve(signal, np.ones(window_size) / window_size, mode="valid")
    variance = np.array(
        [
            np.var(signal[i : i + window_size])
            for i in range(len(signal) - window_size + 1)
        ]
    )

    return moving_avg, variance


def get_maxima_in_moving_interval(
    function_values: npt.NDArray, interval_size: int, step_size: int, filter_value: int
) -> npt.NDArray:
    """
    Slide an interval along the function, filtering out points below a threshold.

    Points smaller than filter_value times the maximum within the interval are removed.

    Parameters
    ----------
    function_values : NDArray
        1D array of function values.
    interval_size : int
        Size of the interval (number of points).
    step_size : int
        Step size for moving the interval.

    Returns
    -------
    NDArray
        Filtered array where points outside the threshold are set to NaN.
    """
    # Convert input to a NumPy array
    function_values = np.asarray(function_values)

    # Initialize an array to store the result
    filtered_indices = []

    # Slide the interval along the function
    for start in range(0, len(function_values) - interval_size + 1, step_size):
        # Define the end of the interval
        end = start + interval_size

        # Extract the interval
        interval = function_values[start:end]

        # Find the maximum value in the interval
        max_value = np.max(interval)

        # Apply the filter
        threshold = filter_value * max_value

        indices = np.where(interval >= threshold)
        filtered_indices += list(start + indices[0])

    return np.array(list(set(filtered_indices)))


def get_pearson_correlation_coefficient(x: npt.NDArray, y: npt.NDArray) -> np.floating:
    """
    TODO.

    Parameters
    ----------
    x : npt.NDArray
        First array of data points.
    y : npt.NDArray
        Second array of data points.

    Returns
    -------
    floating
        Pearson correlation coefficient between x and y.
    """
    mean_x = np.mean(x)
    mean_y = np.mean(y)
    std_x = np.std(x)
    std_y = np.std(y)

    return np.mean((x - mean_x) * (y - mean_y)) / std_x / std_y


def get_t_test(x: npt.NDArray, y: npt.NDArray) -> np.floating:
    """
    Calculate the t-test statistic for two sets of data.

    Parameters
    ----------
    x : NDArray
        First array of data points.
    y : NDArray
        Second array of data points.

    Returns
    -------
    floating
        t-test statistic for the two sets of data.
    """
    r = get_pearson_correlation_coefficient(x, y)
    n = len(x)

    return np.abs(r) * np.sqrt((n - 2) / (1 - r**2))


def probability_density(t: float, n: int) -> float:
    """
    Probability density function for the t-distribution.

    Parameters
    ----------
    t : float
        The t-value for which the probability density is calculated.
    n : int
        The number of data points in the sample.

    Returns
    -------
    float
        The probability density at the given t-value.
    """
    degrees_of_freedom = n - 2

    return (
        scipy.special.gamma((degrees_of_freedom + 1.0) / 2.0)
        / (
            np.sqrt(np.pi * degrees_of_freedom)
            * scipy.special.gamma(degrees_of_freedom / 2.0)
        )
        * (1 + t**2 / degrees_of_freedom) ** (-(degrees_of_freedom + 1.0) / 2.0)
    )


def get_significance(x: npt.NDArray, t: float) -> float:
    """
    Calculate the significance of a t-test statistic.

    Parameters
    ----------
    x : npt.NDArray
        Array of data points.
    t : float
        t-test statistic value.

    Returns
    -------
    float
        The significance level of the t-test statistic.
    """
    n = len(x)

    return scipy.integrate.quad(probability_density, -np.inf, t, args=(n))[0]


def squared_exponential_kernel(
    x1_vec: npt.NDArray, x2_vec: npt.NDArray, tau: float
) -> npt.NDArray[np.float64]:
    """
    Return a simple squared exponential kernel to determine a similarity measure.

    Parameters
    ----------
    x1_vec : NDArray
        Descriptor for first set of data of shape
        [data points, descriptor dimensions].
    x2_vec : NDArray
        Descriptor for second set of data of shape
        [data points, descriptor dimensions].
    tau : float
        Correlation length for the descriptor.

    Returns
    -------
    NDArray[float64]
        Matrix contianing pairwise kernel values.
    """
    # Ensure inputs are at least 2D (n_samples, n_features)
    x1 = np.atleast_2d(x1_vec)
    x2 = np.atleast_2d(x2_vec)

    # If they are 1D row-vectors, convert to column vectors
    if x1.shape[0] == 1 or x1.shape[1] == 1:
        x1 = x1.reshape(-1, 1)
    if x2.shape[0] == 1 or x2.shape[1] == 1:
        x2 = x2.reshape(-1, 1)

    # Compute squared distances (broadcasting works for both 1D and 2D now)
    diff = x1[:, np.newaxis, :] - x2[np.newaxis, :, :]
    sq_dist = np.sum(diff**2, axis=2)

    # Apply the RBF formula
    return np.exp(-0.5 * sq_dist / tau**2)


class GPR:
    """
    A simple Gaussian Process Regression (GPR) model for interpolation and smoothing.

    Parameters
    ----------
    x : NDArray
        Descriptor of shape [data points, descriptor dimensions].
    y : NDArray
        Data to be learned.
    tau : float
        Correlation length for the descriptor.
    sigma : float
        Uncertainty of the input data.
    """

    def __init__(self, x: npt.NDArray, y: npt.NDArray, tau: float, sigma: float):
        K1 = squared_exponential_kernel(x, x, tau)

        self.K1_inv = np.linalg.inv(K1 + np.eye(len(x)) * sigma)
        self.x = x
        self.y = y
        self.tau = tau
        self.sigma = sigma

    def __call__(self, x: npt.NDArray) -> npt.NDArray:
        return self.predict(x)

    def predict(self, x: npt.NDArray) -> npt.NDArray:
        """
        TODO.

        Parameters
        ----------
        x : NDArray
            TODO

        Returns
        -------
        NDArray
            TODO
        """
        K2 = squared_exponential_kernel(x, self.x, self.tau)
        y_test0 = K2.dot(self.K1_inv)

        return y_test0.dot(self.y)

1	from copy import deepcopy	1✔
2	from typing import Any	1✔
3
4	import numpy as np	1✔
5	import numpy.typing as npt	1✔
6	import scipy	1✔
7
8
9	def get_rotation_matrix(	1✔
10	vec_start: npt.NDArray[np.float64], vec_end: npt.NDArray[np.float64]
11	) -> npt.NDArray[np.float64]:
12	"""
13	Calculate the rotation matrix to align two unit vectors.
14
15	Given a two (unit) vectors, vec_start and vec_end, this function calculates the
16	rotation matrix U, so that U * vec_start = vec_end.
17
18	U the is rotation matrix that rotates vec_start to point in the direction
19	of vec_end.
20
21	https://math.stackexchange.com/questions/180418/calculate-rotation-matrix-to-align-vector-a-to-vector-b-in-3d/897677
22
23	Parameters
24	----------
25	vec_start, vec_end : NDArray[float64]
26	Two vectors that should be aligned. Both vectors must have a l2-norm of 1.
27
28	Returns
29	-------
30	NDArray[float64]
31	The rotation matrix U as npt.NDArray with shape (3,3)
32	"""
33	if not np.isclose(np.linalg.norm(vec_start), 1) and not np.isclose(	1✔
34	np.linalg.norm(vec_end), 1
35	):
36	raise ValueError("`vec_start` and `vec_end` args must be unit vectors")	×
37
38	v = np.cross(vec_start, vec_end)	1✔
39	c = np.dot(vec_start, vec_end)	1✔
40	v_x = np.array([[0, -v[2], v[1]], [v[2], 0, -v[0]], [-v[1], v[0], 0]])	1✔
41	return np.eye(3) + v_x + v_x.dot(v_x) / (1 + c)	1✔
42
43
44	def get_rotation_matrix_around_axis(axis: npt.NDArray, phi: float) -> npt.NDArray:	1✔
45	"""
46	Generate a rotation matrix around a given vector.
47
48	Parameters
49	----------
50	axis : NDArray
51	Axis around which the rotation is done.
52	phi : float
53	Angle of rotation around axis in radiants.
54
55	Returns
56	-------
57	NDArray
58	Rotation matrix
59	"""
60	axis_vec = np.array(axis, dtype=np.float64)	1✔
61	axis_vec /= np.linalg.norm(axis_vec)	1✔
62
63	eye = np.eye(3, dtype=np.float64)	1✔
64	ddt = np.outer(axis_vec, axis_vec)	1✔
65	skew = np.array(	1✔
66	[
67	[0, axis_vec[2], -axis_vec[1]],
68	[-axis_vec[2], 0, axis_vec[0]],
69	[axis_vec[1], -axis_vec[0], 0],
70	],
71	dtype=np.float64,
72	)
73
74	return ddt + np.cos(phi) * (eye - ddt) + np.sin(phi) * skew	1✔
75
76
77	def get_rotation_matrix_around_z_axis(phi: float) -> npt.NDArray:	1✔
78	"""
79	Generate a rotation matrix around the z axis.
80
81	Parameters
82	----------
83	phi : float
84	Angle of rotation around axis in radiants.
85
86	Returns
87	-------
88	NDArray
89	Rotation matrix
90	"""
91	return get_rotation_matrix_around_axis(np.array([0.0, 0.0, 1.0]), phi)	×
92
93
94	def get_mirror_matrix(normal_vector: npt.NDArray) -> npt.NDArray:	1✔
95	"""
96	Generate a transformation matrix for mirroring through plane given by the normal.
97
98	Parameters
99	----------
100	normal_vector : NDArray
101	Normal vector of the mirror plane.
102
103	Returns
104	-------
105	NDArray
106	Mirror matrix
107	"""
108	n_vec = normal_vector / np.linalg.norm(normal_vector)	1✔
109	eps = np.finfo(np.float64).eps	1✔
110	a = n_vec[0]	1✔
111	b = n_vec[1]	1✔
112	c = n_vec[2]	1✔
113	M = np.array(	1✔
114	[
115	[1 - 2 * a*2, -2 a * b, -2 * a * c],
116	[-2 * a * b, 1 - 2 * b*2, -2 b * c],
117	[-2 * a * c, -2 * b * c, 1 - 2 * c**2],
118	]
119	)
120	M[np.abs(M) < eps * 10] = 0	1✔
121	return M	1✔
122
123
124	def get_angle_between_vectors(	1✔
125	vector_1: npt.NDArray, vector_2: npt.NDArray
126	) -> npt.NDArray:
127	"""
128	Determine angle between two vectors.
129
130	Parameters
131	----------
132	vector_1 : NDArray
133	vector_2 : NDArray
134
135	Returns
136	-------
137	float
138	Angle in radiants.
139	"""
140	return (	1✔
141	np.dot(vector_1, vector_2) / np.linalg.norm(vector_1) / np.linalg.norm(vector_2)
142	)
143
144
145	def get_fractional_coords(	1✔
146	cartesian_coords: npt.NDArray, lattice_vectors: npt.NDArray
147	) -> npt.NDArray:
148	"""
149	Transform cartesian coordinates into fractional coordinates.
150
151	Parameters
152	----------
153	cartesian_coords: NDArray
154	Cartesian coordinates of atoms (can be Nx2 or Nx3)
155	lattice_vectors: [N_dim x N_dim] numpy array:
156	Matrix of lattice vectors: Each ROW corresponds to one lattice vector!
157
158	Returns
159	-------
160	fractional_coords: [N x N_dim] numpy array
161	Fractional coordinates of atoms
162	"""
163	fractional_coords = np.linalg.solve(lattice_vectors.T, cartesian_coords.T)	1✔
164	return fractional_coords.T	1✔
165
166
167	def get_cartesian_coords(	1✔
168	frac_coords: npt.NDArray, lattice_vectors: npt.NDArray
169	) -> npt.NDArray:
170	"""
171	Transform fractional coordinates into cartesian coordinates.
172
173	Parameters
174	----------
175	frac_coords: NDArray
176	Fractional coordinates of atoms (can be Nx2 or Nx3)
177	lattice_vectors: NDArray
178	Matrix of lattice vectors: Each ROW corresponds to one lattice vector!
179
180	Returns
181	-------
182	NDArray
183	Cartesian coordinates of atoms
184
185	"""
186	return np.dot(frac_coords, lattice_vectors)	1✔
187
188
189	def get_triple_product(a: npt.NDArray, b: npt.NDArray, c: npt.NDArray) -> npt.NDArray:	1✔
190	"""
191	Get the triple product (DE: Spatprodukt): a*(bxc).
192
193	Parameters
194	----------
195	a: NDArray
196	TODO
197	b: NDArray
198	TODO
199	c: NDArray
200	TODO
201
202	Returns
203	-------
204	NDarray
205	Triple product of each input vector
206	"""
207	if len(a) != 3 or len(b) != 3 or len(c) != 3:	1✔
208	raise ValueError("Each vector must be of length 3")	×
209
210	return np.dot(np.cross(a, b), c)	1✔
211
212
213	def smooth_function(y: npt.NDArray, box_pts: int) -> npt.NDArray[np.floating[Any]]:	1✔
214	"""
215	Smooths a function using convolution.
216
217	Parameters
218	----------
219	y : NDArray
220	TODO
221	box_pts : int
222	TODO
223
224	Returns
225	-------
226	y_smooth : NDArray[floating[Any]]
227	TODO
228	"""
229	box = np.ones(box_pts) / box_pts	1✔
230	return np.convolve(y, box, mode="same")	1✔
231
232
233	def get_cross_correlation_function(	1✔
234	signal_0: npt.NDArray,
235	signal_1: npt.NDArray,
236	detrend: bool = False,
237	) -> npt.NDArray:
238	"""
239	Calculate the autocorrelation function for a given signal.
240
241	Parameters
242	----------
243	signal_0 : NDArray
244	First siganl for which the correlation function should be calculated.
245	signal_1 : NDArray
246	Second siganl for which the correlation function should be calculated.
247
248	Returns
249	-------
250	NDArray
251	Autocorrelation function from 0 to max_lag.
252	"""
253	if detrend:	1✔
254	signal_0 = scipy.signal.detrend(signal_0)	×
255	signal_1 = scipy.signal.detrend(signal_1)	×
256
257	cross_correlation = np.correlate(signal_0, signal_1, mode="full")	1✔
258	cross_correlation = cross_correlation[cross_correlation.size // 2 :]	1✔
259
260	# normalize by number of overlapping data points
261	cross_correlation /= np.arange(cross_correlation.size, 0, -1)	1✔
262	cutoff = int(cross_correlation.size * 0.75)	1✔
263	return cross_correlation[:cutoff]	1✔
264
265
266	def get_autocorrelation_function_manual_lag(	1✔
267	signal: npt.NDArray, max_lag: int
268	) -> npt.NDArray:
269	"""
270	Alternative method for autocorrelation of a signal using numpy.corrcoef.
271
272	Allows the lag to be set manually.
273
274	Parameters
275	----------
276	signal : NDArray
277	Siganl for which the autocorrelation function should be calculated.
278	max_lag : int \| None
279	Autocorrelation will be calculated for a range of 0 to max_lag,
280	where max_lag is the largest lag for the calculation of the
281	autocorrelation function
282
283	Returns
284	-------
285	autocorrelation : npt.NDArray
286	Autocorrelation function from 0 to max_lag.
287	"""
288	lag = npt.NDArray(range(max_lag))	×
289	autocorrelation = np.array([np.nan] * max_lag)	×
290
291	for i in lag:	×
292	corr = 1.0 if i == 0 else np.corrcoef(signal[i], signal[:-i])[0][1]	×
293	autocorrelation[i] = corr	×
294
295	return autocorrelation	×
296
297
298	def get_fourier_transform(	1✔
299	signal: npt.NDArray, time_step: float
300	) -> tuple[npt.NDArray, npt.NDArray]:
301	"""
302	Calculate the fourier transform of a given siganl.
303
304	Parameters
305	----------
306	signal : NDArray
307	Siganl for which the autocorrelation function should be calculated.
308	time_step : float
309	Time step of the signal in seconds.
310
311	Returns
312	-------
313	tuple[NDArray, NDArray]
314	Frequencs and absolute values of the fourier transform.
315
316	"""
317	f = scipy.fft.fftfreq(signal.size, d=time_step)	1✔
318	y = scipy.fft.fft(signal)	1✔
319
320	L = f >= 0	1✔
321
322	return f[L], y[L]	1✔
323
324
325	def lorentzian(	1✔
326	x: float \| npt.NDArray[np.float64], a: float, b: float, c: float
327	) -> float \| npt.NDArray[np.float64]:
328	"""
329	Return a Lorentzian function.
330
331	Parameters
332	----------
333	x : float \| NDArray[float64]
334	Argument x of f(x) --> y.
335	a : float
336	Maximum of Lorentzian.
337	b : float
338	Width of Lorentzian.
339	c : float
340	Magnitude of Lorentzian.
341
342	Returns
343	-------
344	f : float \| NDArray[float64]
345	Outupt of a Lorentzian function.
346	"""
347	return c / (1.0 + ((x - a) / (b / 2.0)) ** 2)	1✔
348
349
350	def exponential(	1✔
351	x: float \| npt.NDArray[np.float64], a: float, b: float
352	) -> float \| npt.NDArray[np.float64]:
353	"""
354	Return an exponential function.
355
356	Parameters
357	----------
358	x : float \| NDArray[float64]
359	Argument x of f(x) --> y
360	a : float
361	TODO
362	b : float
363	TODO
364	"""
365	return a * np.exp(x * b)	×
366
367
368	def double_exponential(	1✔
369	x: float \| npt.NDArray, a: float, b: float, c: float
370	) -> float \| npt.NDArray:
371	"""TODO."""
372	return a * (np.exp(x * b) + np.exp(x * c))	×
373
374
375	def gaussian_window(N: int, std: float = 0.4) -> npt.NDArray[np.float64]:	1✔
376	"""
377	Generate a Gaussian window.
378
379	Parameters
380	----------
381	N : int
382	Number of points in the window.
383	std : float
384	Standard deviation of the Gaussian window, normalized
385	such that the maximum value occurs at the center of the window.
386
387	Returns
388	-------
389	window : NDarray[float64]
390	Gaussian window of length N.
391
392	"""
393	n = np.linspace(-1, 1, N)	1✔
394	return np.exp(-0.5 * (n / std) ** 2)	1✔
395
396
397	def apply_gaussian_window(	1✔
398	data: npt.NDArray[np.float64], std: float = 0.4
399	) -> npt.NDArray[np.float64]:
400	"""
401	Apply a Gaussian window to an array.
402
403	Parameters
404	----------
405	data : NDarray[float64]
406	Input data array to be windowed.
407	std : float
408	Standard deviation of the Gaussian window.
409
410	Returns
411	-------
412	windowed_data : NDArray[float64]
413	Windowed data array.
414	"""
415	N = len(data)	1✔
416	window = gaussian_window(N, std)	1✔
417	return data * window	1✔
418
419
420	def hann_window(N: int) -> npt.NDArray[np.float64]:	1✔
421	"""
422	Generate a Hann window.
423
424	Parameters
425	----------
426	N : int
427	Number of points in the window.
428
429	Returns
430	-------
431	NDArray
432	Hann window of length N.
433	"""
434	return 0.5 * (1 - np.cos(2 * np.pi * np.arange(N) / (N - 1)))	×
435
436
437	def apply_window(data: npt.NDArray[np.float64]) -> npt.NDArray[np.float64]:	1✔
438	"""
439	Apply a Hann window to an array.
440
441	Parameters
442	----------
443	data : NDarray[float64]
444	Input data array to be windowed.
445
446	Returns
447	-------
448	NDarray[float64]
449	Windowed data array.
450	"""
451	N = len(data)	×
452	window = hann_window(N)	×
453	return data * window	×
454
455
456	def norm_matrix_by_dagonal(matrix: npt.NDArray) -> npt.NDArray:	1✔
457	"""
458	Norms a matrix such that the diagonal becomes 1.
459
460	\| a_11 a_12 a_13 \| \| 1 a'_12 a'_13 \|
461	\| a_21 a_22 a_23 \| --> \| a'_21 1 a'_23 \|
462	\| a_31 a_32 a_33 \| \| a'_31 a'_32 1 \|
463
464	Parameters
465	----------
466	matrix : NDArray
467	Matrix that should be normed.
468
469	Returns
470	-------
471	matrix : NDArray
472	Normed matrix.
473	"""
474	diagonal = np.array(np.diagonal(matrix))	1✔
475	L = diagonal == 0.0	1✔
476	diagonal[L] = 1.0	1✔
477
478	new_matrix = deepcopy(matrix)	1✔
479	new_matrix /= np.sqrt(	1✔
480	np.tile(diagonal, (matrix.shape[1], 1)).T
481	* np.tile(diagonal, (matrix.shape[0], 1))
482	)
483
484	return new_matrix	1✔
485
486
487	def mae(delta: npt.NDArray) -> np.floating:	1✔
488	"""
489	Calculate the mean absolute error from a list of value differnces.
490
491	Parameters
492	----------
493	delta : NDArray
494	Array containing differences
495
496	Returns
497	-------
498	float
499	mean absolute error
500	"""
501	return np.mean(np.abs(delta))	1✔
502
503
504	def rel_mae(	1✔
505	delta: npt.NDArray, target_val: npt.NDArray
506	) -> np.floating \| np.float64 \| float:
507	"""
508	Compute relative MAE from value differences and corresponding target values.
509
510	Parameters
511	----------
512	delta : NDArray
513	Array containing differences
514	target_val : NDArray
515	Array of target values against which the difference should be compared
516
517	Returns
518	-------
519	float
520	relative mean absolute error
521	"""
522	target_norm = np.mean(np.abs(target_val))	1✔
523	return np.mean(np.abs(delta)).item() / (target_norm + 1e-9)	1✔
524
525
526	def rmse(delta: npt.NDArray) -> float:	1✔
527	"""
528	Calculate the root mean sqare error from a list of value differnces.
529
530	Parameters
531	----------
532	delta : NDArray
533	Array containing differences
534
535	Returns
536	-------
537	float
538	root mean square error
539	"""
540	return np.sqrt(np.mean(np.square(delta)))	×
541
542
543	def rel_rmse(delta: npt.NDArray, target_val: npt.NDArray) -> float:	1✔
544	"""
545	Calculate the relative root mean sqare error from a list of value differences.
546
547	Parameters
548	----------
549	delta : NDArray
550	Array containing differences
551	target_val : NDArray
552	Array of target values against which the difference should be compared
553
554	Returns
555	-------
556	float
557	relative root mean sqare error
558
559	"""
560	target_norm = np.sqrt(np.mean(np.square(target_val)))	×
561	return np.sqrt(np.mean(np.square(delta))) / (target_norm + 1e-9)	×
562
563
564	def get_moving_average(	1✔
565	signal: npt.NDArray[np.float64], window_size: int
566	) -> tuple[npt.NDArray[np.floating], npt.NDArray[np.floating]]:
567	"""
568	Cacluated the moving average and the variance around the moving average.
569
570	Parameters
571	----------
572	signal : NDArray[float64]
573	Signal for which the moving average should be calculated.
574	window_size : int
575	Window size for the mocing average.
576
577	Returns
578	-------
579	npt.NDArray[floating]
580	Moving average.
581	npt.NDArray[floating]
582	Variance around the moving average.
583	"""
584	moving_avg = np.convolve(signal, np.ones(window_size) / window_size, mode="valid")	×
585	variance = np.array(	×
586	[
587	np.var(signal[i : i + window_size])
588	for i in range(len(signal) - window_size + 1)
589	]
590	)
591
592	return moving_avg, variance	×
593
594
595	def get_maxima_in_moving_interval(	1✔
596	function_values: npt.NDArray, interval_size: int, step_size: int, filter_value: int
597	) -> npt.NDArray:
598	"""
599	Slide an interval along the function, filtering out points below a threshold.
600
601	Points smaller than filter_value times the maximum within the interval are removed.
602
603	Parameters
604	----------
605	function_values : NDArray
606	1D array of function values.
607	interval_size : int
608	Size of the interval (number of points).
609	step_size : int
610	Step size for moving the interval.
611
612	Returns
613	-------
614	NDArray
615	Filtered array where points outside the threshold are set to NaN.
616	"""
617	# Convert input to a NumPy array
618	function_values = np.asarray(function_values)	×
619
620	# Initialize an array to store the result
621	filtered_indices = []	×
622
623	# Slide the interval along the function
624	for start in range(0, len(function_values) - interval_size + 1, step_size):	×
625	# Define the end of the interval
626	end = start + interval_size	×
627
628	# Extract the interval
629	interval = function_values[start:end]	×
630
631	# Find the maximum value in the interval
632	max_value = np.max(interval)	×
633
634	# Apply the filter
635	threshold = filter_value * max_value	×
636
637	indices = np.where(interval >= threshold)	×
638	filtered_indices += list(start + indices[0])	×
639
640	return np.array(list(set(filtered_indices)))	×
641
642
643	def get_pearson_correlation_coefficient(x: npt.NDArray, y: npt.NDArray) -> np.floating:	1✔
644	"""
645	TODO.
646
647	Parameters
648	----------
649	x : npt.NDArray
650	First array of data points.
651	y : npt.NDArray
652	Second array of data points.
653
654	Returns
655	-------
656	floating
657	Pearson correlation coefficient between x and y.
658	"""
659	mean_x = np.mean(x)	×
660	mean_y = np.mean(y)	×
661	std_x = np.std(x)	×
662	std_y = np.std(y)	×
663
664	return np.mean((x - mean_x) * (y - mean_y)) / std_x / std_y	×
665
666
667	def get_t_test(x: npt.NDArray, y: npt.NDArray) -> np.floating:	1✔
668	"""
669	Calculate the t-test statistic for two sets of data.
670
671	Parameters
672	----------
673	x : NDArray
674	First array of data points.
675	y : NDArray
676	Second array of data points.
677
678	Returns
679	-------
680	floating
681	t-test statistic for the two sets of data.
682	"""
683	r = get_pearson_correlation_coefficient(x, y)	×
684	n = len(x)	×
685
686	return np.abs(r) * np.sqrt((n - 2) / (1 - r**2))	×
687
688
689	def probability_density(t: float, n: int) -> float:	1✔
690	"""
691	Probability density function for the t-distribution.
692
693	Parameters
694	----------
695	t : float
696	The t-value for which the probability density is calculated.
697	n : int
698	The number of data points in the sample.
699
700	Returns
701	-------
702	float
703	The probability density at the given t-value.
704	"""
705	degrees_of_freedom = n - 2	×
706
707	return (	×
708	scipy.special.gamma((degrees_of_freedom + 1.0) / 2.0)
709	/ (
710	np.sqrt(np.pi * degrees_of_freedom)
711	* scipy.special.gamma(degrees_of_freedom / 2.0)
712	)
713	* (1 + t2 / degrees_of_freedom) (-(degrees_of_freedom + 1.0) / 2.0)
714	)
715
716
717	def get_significance(x: npt.NDArray, t: float) -> float:	1✔
718	"""
719	Calculate the significance of a t-test statistic.
720
721	Parameters
722	----------
723	x : npt.NDArray
724	Array of data points.
725	t : float
726	t-test statistic value.
727
728	Returns
729	-------
730	float
731	The significance level of the t-test statistic.
732	"""
733	n = len(x)	×
734
735	return scipy.integrate.quad(probability_density, -np.inf, t, args=(n))[0]	×
736
737
738	def squared_exponential_kernel(	1✔
739	x1_vec: npt.NDArray, x2_vec: npt.NDArray, tau: float
740	) -> npt.NDArray[np.float64]:
741	"""
742	Return a simple squared exponential kernel to determine a similarity measure.
743
744	Parameters
745	----------
746	x1_vec : NDArray
747	Descriptor for first set of data of shape
748	[data points, descriptor dimensions].
749	x2_vec : NDArray
750	Descriptor for second set of data of shape
751	[data points, descriptor dimensions].
752	tau : float
753	Correlation length for the descriptor.
754
755	Returns
756	-------
757	NDArray[float64]
758	Matrix contianing pairwise kernel values.
759	"""
760	# Ensure inputs are at least 2D (n_samples, n_features)
761	x1 = np.atleast_2d(x1_vec)	×
762	x2 = np.atleast_2d(x2_vec)	×
763
764	# If they are 1D row-vectors, convert to column vectors
765	if x1.shape[0] == 1 or x1.shape[1] == 1:	×
766	x1 = x1.reshape(-1, 1)	×
767	if x2.shape[0] == 1 or x2.shape[1] == 1:	×
768	x2 = x2.reshape(-1, 1)	×
769
770	# Compute squared distances (broadcasting works for both 1D and 2D now)
771	diff = x1[:, np.newaxis, :] - x2[np.newaxis, :, :]	×
772	sq_dist = np.sum(diff**2, axis=2)	×
773
774	# Apply the RBF formula
775	return np.exp(-0.5 * sq_dist / tau**2)	×
776
777
778	class GPR:	1✔
779	"""
780	A simple Gaussian Process Regression (GPR) model for interpolation and smoothing.
781
782	Parameters
783	----------
784	x : NDArray
785	Descriptor of shape [data points, descriptor dimensions].
786	y : NDArray
787	Data to be learned.
788	tau : float
789	Correlation length for the descriptor.
790	sigma : float
791	Uncertainty of the input data.
792	"""
793
794	def __init__(self, x: npt.NDArray, y: npt.NDArray, tau: float, sigma: float):	1✔
795	K1 = squared_exponential_kernel(x, x, tau)	×
796
797	self.K1_inv = np.linalg.inv(K1 + np.eye(len(x)) * sigma)	×
798	self.x = x	×
799	self.y = y	×
800	self.tau = tau	×
801	self.sigma = sigma	×
802
803	def __call__(self, x: npt.NDArray) -> npt.NDArray:	1✔
804	return self.predict(x)	×
805
806	def predict(self, x: npt.NDArray) -> npt.NDArray:	1✔
807	"""
808	TODO.
809
810	Parameters
811	----------
812	x : NDArray
813	TODO
814
815	Returns
816	-------
817	NDArray
818	TODO
819	"""
820	K2 = squared_exponential_kernel(x, self.x, self.tau)	×
821	y_test0 = K2.dot(self.K1_inv)	×
822
823	return y_test0.dot(self.y)	×

maurergroup / dfttoolkit / 15472938725

Source File Press 'n' to go to next uncovered line, 'b' for previous

Source File
Press 'n' to go to next uncovered line, 'b' for previous