import numpy as np

from preprocessing.constants import CA_C_DIST, N_CA_DIST, N_CA_C_ANGLE


def rotation_matrix(angle, axis):
    """
    Args:
        angle: (n,)
        axis: 0=x, 1=y, 2=z
    Returns:
        (n, 3, 3)
    """
    n = len(angle)
    R = np.eye(3)[None, :, :].repeat(n, axis=0)

    axis = 2 - axis
    start = axis // 2
    step = axis % 2 + 1
    s = slice(start, start + step + 1, step)

    R[:, s, s] = np.array(
        [[np.cos(angle), (-1) ** (axis + 1) * np.sin(angle)],
         [(-1) ** axis * np.sin(angle), np.cos(angle)]]
    ).transpose(2, 0, 1)
    return R


def get_bb_transform(n_xyz, ca_xyz, c_xyz):
    """
    Compute translation and rotation of the canoncical backbone frame (triangle N-Ca-C) from a position with
    Ca at the origin, N on the x-axis and C in the xy-plane to the global position of the backbone frame

    Args:
        n_xyz: (n, 3)
        ca_xyz: (n, 3)
        c_xyz: (n, 3)

    Returns:
        quaternion represented as array of shape (n, 4)
        translation vector which is an array of shape (n, 3)
    """

    translation = ca_xyz
    n_xyz = n_xyz - translation
    c_xyz = c_xyz - translation

    # Find rotation matrix that aligns the coordinate systems
    #    rotate around y-axis to move N into the xy-plane
    theta_y = np.arctan2(n_xyz[:, 2], -n_xyz[:, 0])
    Ry = rotation_matrix(theta_y, 1)
    n_xyz = np.einsum('noi,ni->no', Ry.transpose(0, 2, 1), n_xyz)

    #    rotate around z-axis to move N onto the x-axis
    theta_z = np.arctan2(n_xyz[:, 1], n_xyz[:, 0])
    Rz = rotation_matrix(theta_z, 2)
    # n_xyz = np.einsum('noi,ni->no', Rz.transpose(0, 2, 1), n_xyz)

    #    rotate around x-axis to move C into the xy-plane
    c_xyz = np.einsum('noj,nji,ni->no', Rz.transpose(0, 2, 1),
                      Ry.transpose(0, 2, 1), c_xyz)
    theta_x = np.arctan2(c_xyz[:, 2], c_xyz[:, 1])
    Rx = rotation_matrix(theta_x, 0)

    # Final rotation matrix
    R = np.einsum('nok,nkj,nji->noi', Ry, Rz, Rx)

    # Convert to quaternion
    # q = w + i*u_x + j*u_y + k * u_z
    quaternion = rotation_matrix_to_quaternion(R)

    return quaternion, translation


def get_bb_coords_from_transform(ca_coords, quaternion):
    """
    Args:
        ca_coords: (n, 3)
        quaternion: (n, 4)
    Returns:
        backbone coords (n*3, 3), order is [N, CA, C]
        backbone atom types as a list of length n*3
    """
    R = quaternion_to_rotation_matrix(quaternion)
    bb_coords = np.tile(np.array(
        [[N_CA_DIST, 0, 0],
         [0, 0, 0],
         [CA_C_DIST * np.cos(N_CA_C_ANGLE), CA_C_DIST * np.sin(N_CA_C_ANGLE), 0]]),
        [len(ca_coords), 1])
    bb_coords = np.einsum('noi,ni->no', R.repeat(3, axis=0), bb_coords) + ca_coords.repeat(3, axis=0)
    bb_atom_types = [t for _ in range(len(ca_coords)) for t in ['N', 'C', 'C']]

    return bb_coords, bb_atom_types


def quaternion_to_rotation_matrix(q):
    """
    x_rot = R x

    Args:
        q: (n, 4)
    Returns:
        R: (n, 3, 3)
    """
    # Normalize
    q = q / (q ** 2).sum(1, keepdims=True) ** 0.5

    # https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
    w, x, y, z = q[:, 0], q[:, 1], q[:, 2], q[:, 3]
    R = np.stack([
        np.stack([1 - 2 * y ** 2 - 2 * z ** 2, 2 * x * y - 2 * z * w,
                  2 * x * z + 2 * y * w], axis=1),
        np.stack([2 * x * y + 2 * z * w, 1 - 2 * x ** 2 - 2 * z ** 2,
                  2 * y * z - 2 * x * w], axis=1),
        np.stack([2 * x * z - 2 * y * w, 2 * y * z + 2 * x * w,
                  1 - 2 * x ** 2 - 2 * y ** 2], axis=1)
    ], axis=1)

    return R


def rotation_matrix_to_quaternion(R):
    """
    https://en.wikipedia.org/wiki/Rotation_matrix#Quaternion
    Args:
        R: (n, 3, 3)
    Returns:
        q: (n, 4)
    """

    t = R[:, 0, 0] + R[:, 1, 1] + R[:, 2, 2]
    r = np.sqrt(1 + t)
    w = 0.5 * r
    x = np.sign(R[:, 2, 1] - R[:, 1, 2]) * np.abs(
        0.5 * np.sqrt(1 + R[:, 0, 0] - R[:, 1, 1] - R[:, 2, 2]))
    y = np.sign(R[:, 0, 2] - R[:, 2, 0]) * np.abs(
        0.5 * np.sqrt(1 - R[:, 0, 0] + R[:, 1, 1] - R[:, 2, 2]))
    z = np.sign(R[:, 1, 0] - R[:, 0, 1]) * np.abs(
        0.5 * np.sqrt(1 - R[:, 0, 0] - R[:, 1, 1] + R[:, 2, 2]))

    return np.stack((w, x, y, z), axis=1)
