14713210967

from copy import deepcopy
from typing import Union

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


def get_rotation_matrix(
    vec_start: npt.NDArray, vec_end: npt.NDArray
) -> npt.NDArray:
    """
    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 : npt.NDArray[np.float64]
        Two vectors that should be aligned. Both vectors must have a l2-norm of 1.

    Returns:
    --------
    R
        The rotation matrix U as npt.NDArray with shape (3,3)
    """
    assert np.isclose(np.linalg.norm(vec_start), 1) and np.isclose(
        np.linalg.norm(vec_end), 1
    ), "vec_start and vec_end 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]])
    R = np.eye(3) + v_x + v_x.dot(v_x) / (1 + c)

    return R


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

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

    Returns
    -------
    R : npt.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,
    )

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


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

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

    Returns
    -------
    npt.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:
    """
    Generates a transformation matrix for mirroring through plane given by the
    normal vector.

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

    Returns
    -------
    M : npt.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:
    """
    Determines angle between two vectors.

    Parameters
    ----------
    vector_1 : npt.NDArray
    vector_2 : npt.NDArray

    Returns
    -------
    angle : float
        Angle in radiants.

    """
    angle = (
        np.dot(vector_1, vector_2)
        / np.linalg.norm(vector_1)
        / np.linalg.norm(vector_2)
    )
    return angle


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

    Parameters
    ----------
    cartesian_coords: [N x N_dim] numpy array
        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: [N x N_dim] numpy array
        Fractional 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
    -------
    cartesian_coords: [N x N_dim] numpy array
        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:
    """
    Returns the triple product (DE: Spatprodukt): a*(bxc)

    """
    assert len(a) == 3
    assert len(b) == 3
    assert len(c) == 3
    return np.dot(np.cross(a, b), c)


def smooth_function(y: npt.NDArray, box_pts: int):
    """
    Smooths a function using convolution.

    Parameters
    ----------
    y : TYPE
        DESCRIPTION.
    box_pts : TYPE
        DESCRIPTION.

    Returns
    -------
    y_smooth : TYPE
        DESCRIPTION.

    """
    box = np.ones(box_pts) / box_pts
    y_smooth = np.convolve(y, box, mode="same")
    return y_smooth


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 : 1D npt.NDArray
        First siganl for which the correlation function should be calculated.
    signal_1 : 1D npt.NDArray
        Second siganl for which the correlation function should be calculated.

    Returns
    -------
    correlation : npt.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='same')
    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)
    cross_correlation = cross_correlation[:cutoff]

    return cross_correlation


def get_autocorrelation_function_manual_lag(
    signal: npt.NDArray, max_lag: int
) -> npt.NDArray:
    """
    Alternative method to determine the autocorrelation function for a given
    signal that used numpy.corrcoef. This function allows to set the lag
    manually.

    Parameters
    ----------
    signal : 1D npt.NDArray
        Siganl for which the autocorrelation function should be calculated.
    max_lag : Union[None, int]
        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 l in lag:
        if l == 0:
            corr = 1.0
        else:
            corr = np.corrcoef(signal[l:], signal[:-l])[0][1]

        autocorrelation[l] = corr

    return autocorrelation


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

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

    Returns
    -------
    (npt.NDArray, npt.NDArray)
        Frequencs and absolute values of the fourier transform.

    """
    # d = len(signal) * time_step

    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: Union[float, npt.NDArray], a: float, b: float, c: float
) -> Union[float, npt.NDArray]:
    """
    Returns a Lorentzian function.

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

    Returns
    -------
    f : Union[float, npt.NDArray]
        Outupt of a Lorentzian function.

    """
    # f = c / (np.pi * b * (1.0 + ((x - a) / b) ** 2))  # +d
    f = c / (1.0 + ((x - a) / (b / 2.0)) ** 2)  # +d

    return f


def exponential(
    x: Union[float, npt.NDArray], a: float, b: float
) -> Union[float, npt.NDArray]:
    f = a * np.exp(x * b)
    return f


def double_exponential(
    x: Union[float, npt.NDArray], a: float, b: float, c: float
) -> Union[float, npt.NDArray]:
    f = a * (np.exp(x * b) + np.exp(x * c))
    return f


def gaussian_window(N, std=0.4):
    """
    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 : np.array
        Gaussian window of length N.

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


def apply_gaussian_window(data, std=0.4):
    """
    Apply a Gaussian window to an array.

    Parameters
    ----------
    data : np.array
        Input data array to be windowed.
    std : float
        Standard deviation of the Gaussian window.

    Returns
    -------
    windowed_data : np.array
        Windowed data array.

    """
    N = len(data)
    window = gaussian_window(N, std)
    windowed_data = data * window
    return windowed_data


def hann_window(N):
    """
    Generate a Hann window.

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

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


def apply_hann_window(data):
    """
    Apply a Hann window to an array.

    Parameters
    ----------
    data : np.ndarray
        Input data array to be windowed.

    Returns
    -------
    np.ndarray
        Windowed data array.
    """
    N = len(data)
    window = hann_window(N)
    windowed_data = data * window
    return windowed_data


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 : npt.NDArray
        Matrix that should be normed.

    Returns
    -------
    matrix : npt.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: np.ndarray) -> np.floating:
    """
    Calculated the mean absolute error from a list of value differnces.

    Parameters
    ----------
    delta : np.ndarray
        Array containing differences

    Returns
    -------
    float
        mean absolute error

    """
    return np.mean(np.abs(delta))


def rel_mae(delta: np.ndarray, target_val: np.ndarray) -> np.floating:
    """
    Calculated the relative mean absolute error from a list of value differnces,
    given the target values.

    Parameters
    ----------
    delta : np.ndarray
        Array containing differences
    target_val : np.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: np.ndarray) -> float:
    """
    Calculated the root mean sqare error from a list of value differnces.

    Parameters
    ----------
    delta : np.ndarray
        Array containing differences

    Returns
    -------
    float
        root mean square error

    """
    return np.sqrt(np.mean(np.square(delta)))


def rel_rmse(delta: np.ndarray, target_val: np.ndarray) -> float:
    """
    Calculated the relative root mean sqare error from a list of value differnces,
    given the target values.

    Parameters
    ----------
    delta : np.ndarray
        Array containing differences
    target_val : np.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):
    """
    Cacluated the moving average and the variance around the moving average.

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

    Returns
    -------
    moving_avg : TYPE
        Moving average.
    variance : TYPE
        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, interval_size, step_size, filter_value
):
    """
    Moves an interval along the function and filters out all points smaller
    than filter_value time the maximal value in the interval.

    Parameters
    ----------
    function_values : array-like
        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
    -------
        np.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, 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, y):
    r = get_pearson_correlation_coefficient(x, y)

    n = len(x)

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

    return t


def probability_density(t, n):

    degrees_of_freedom = n - 2

    f = (
        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)
    )

    return f


def get_significance(x, t):

    n = len(x)

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

    return Df


# def squared_exponential_kernel(x1_vec_0, x2_vec_0, tau):
#     x1_vec = np.atleast_3d(x1_vec_0)
#     x2_vec = np.atleast_3d(x2_vec_0)

#     x1 = np.tile(x1_vec, len(x2_vec))
#     x2 = np.tile(x2_vec, len(x1_vec))

#     print(x1_vec.shape, x2_vec.shape, x1.shape, x2.shape)

#     x1 = np.moveaxis(x1, 2, 1)
#     x2 = x2.T
#     x2 = np.moveaxis(x2, 2, 1)

#     x_diff = x1 - x2

#     a = np.sum(x_diff**2, axis=2)

#     M = np.exp(-0.5 / tau**2 * a)

#     return M


def squared_exponential_kernel(x1_vec, x2_vec, tau):
    # 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
    K = np.exp(-0.5 * sq_dist / tau**2)

    return K


class GPR:
    def __init__(self, x, y, tau, sigma):
        """


        Parameters
        ----------
        x : np.array
            Descriptor of shape [data points, descriptor dimensions].
        y : TYPE
            DESCRIPTION.
        tau : TYPE
            DESCRIPTION.
        sigma : TYPE
            DESCRIPTION.

        Returns
        -------
        None.

        """
        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):
        return self.predict(x)

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

maurergroup / dfttoolkit / 14713210967

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