Source code for ariane.lib.utils.math

import numpy as np
from scipy import linalg as la, interpolate
from scipy.spatial import distance as dist




def map_intervals(bounds_src, bounds_tar):
    """Return function mapping between two intervals."""
    a_src, b_src = bounds_src
    a_tar, b_tar = bounds_tar

    def f_map(x):
        return a_tar + ((b_tar-a_tar) / (b_src-a_src)) * (x-a_src)
    return f_map





def normalizer(a_src, b_src):
    """Return function mapping an interval to (0, 1)."""
    return map_intervals((a_src, b_src), (0, 1))





def denormalizer(a_src, b_src):
    """Return function mapping (0, 1) to an interval."""
    return map_intervals((0, 1), (a_src, b_src))





def cartesian_product(*s):
    """Return the cartesian product of a set of sets."""
    meshgrid = np.array(np.meshgrid(*s))
    return np.transpose([np.ravel(meshgrid_dim)
                         for meshgrid_dim in meshgrid])





def cholesky_update(L, Aij_new, Aii_new):
    """Update a Cholesky decomposition for new entries."""
    V = la.solve_triangular(L, Aij_new, lower=True).T
    L_ = la.cholesky(Aii_new - V@V.T, lower=True)

    return np.vstack((np.hstack((L, np.zeros_like(Aij_new))),
                      np.hstack((V, L_))))





def compute_corner_pairs(*limits):
    """Return all pairs of corner points of a hypercube."""
    dim = len(limits)
    x_corners = cartesian_product(*limits)
    return x_corners[cartesian_product(*[np.arange(2**dim)]*dim)]





def make_grid(limits, num):
    """Return a grid in a hyperrectangle with given number(s) of points per dimension."""
    limits = np.atleast_2d(limits)
    if type(num) is int:
        num = len(limits)*[num]
    num = np.atleast_1d(num)
    assert len(limits) == len(num)

    locs_per_dim = [np.linspace(limits[i, 0], limits[i, 1], num=num[i])
                    for i in range(len(limits))]
    return cartesian_product(*locs_per_dim)





def rotation_matrix_2d(alpha: float, in_rad=False):
    """Return a 2D matrix rotating a coordinate system by a given angle."""
    if not in_rad:
        alpha = np.deg2rad(alpha)
    c = np.cos(alpha)
    s = np.sin(alpha)
    return np.array([[c, -s], [s, c]])





def normalized_linear_interpolator(x, f, interpolator=interpolate.LinearNDInterpolator,
                                   fill_value=0, x_min=None, x_max=None):
    """Return a linear interpolator interpolating on a normalized (unit) hypercube."""
    x = np.asarray(x)

    assert len(x) == len(f)

    if x_min is None:
        x_min = np.min(x, axis=0)

    if x_max is None:
        x_max = np.max(x, axis=0)

    x_min = np.atleast_1d(x_min)
    x_max = np.atleast_1d(x_max)

    assert np.shape(x)[1] == len(x_min) == len(x_max)
    assert np.alltrue(x_min <= x_max)

    x_normalizers = [normalizer(x_min[i], x_max[i]) for i in range(len(x_min))]

    x_norm = np.transpose([x_normalizers[i](x[:, i]) for i in range(len(x_normalizers))])

    interp = interpolator(x_norm, f, fill_value=fill_value)

    return lambda x: interp([x_normalizers[i](x[i]) for i in range(len(x_normalizers))])[0]