import torch


def rotation_from_sphere_points_torch(x, y):
    """
    Gets the rotation matrix that moves x to y in the geodesic path on the sphere.
    Based on the equations of "Analysis of principal nested spheres", Sung et al. 2012 (appendix)

    Parameters
    ----------
    :param x: point on a sphere
    :param y: point on a sphere

    Returns
    -------
    :return: rotation matrix
    """
    if x.dim() == 1:
        x = x.unsqueeze(-2)
    if y.dim() == 1:
        y = y.unsqueeze(-2)

    dim = x.shape[1]

    # Compute the inner product
    inner_product = torch.mm(x, y.T)
    # Clamp in case any value is not in the interval [-1,1]
    # A small number is added/substracted to the bounds to avoid NaNs during backward computation.
    inner_product = inner_product.clamp(-1. + 1e-15, 1. - 1e-15)

    # Compute intermediate vector
    c_vec = x - y * inner_product
    c_vec = c_vec / torch.norm(c_vec)

    R = torch.eye(dim, dim, dtype=inner_product.dtype) + \
        torch.sin(torch.acos(inner_product)) * (torch.mm(y.T, c_vec) - torch.mm(c_vec.T, y)) + \
        (inner_product - 1.) * (torch.mm(y.T, y) + torch.mm(c_vec.T, c_vec))

    return R