# Copyright (C) 2022, Joao Rodrigues (j.p.g.l.m.rodrigues@gmail.com
#                     Anuj Sharma (anuj.sharma80@gmail.com)
#
# This file is part of the Biopython distribution and governed by your
# choice of the "Biopython License Agreement" or the "BSD 3-Clause License".
# Please see the LICENSE file that should have been included as part of this
# package.

"""Structural alignment using Quaternion Characteristic Polynomial (QCP).

QCPSuperimposer finds the best rotation and translation to put
two point sets on top of each other (minimizing the RMSD). This is
eg. useful to superimpose crystal structures. QCP stands for
Quaternion Characteristic Polynomial, which is used in the algorithm.

Algorithm and original code described in:

Theobald DL.
Rapid calculation of RMSDs using a quaternion-based characteristic polynomial.
Acta Crystallogr A. 2005 Jul;61(Pt 4):478-80. doi: 10.1107/S0108767305015266.
Epub 2005 Jun 23. PMID: 15973002.
"""

import numpy as np

from Bio.PDB.PDBExceptions import PDBException


def qcp(coords1, coords2, natoms):
    """Implement the QCP code in Python.

    Input coordinate arrays must be centered at the origin and have
    shape Nx3.

    Variable names match (as much as possible) the C implementation.
    """
    # Original code has coords1 be the mobile. I think it makes more sense
    # for it to be the reference, so I swapped here.
    G1 = np.trace(np.dot(coords2, coords2.T))
    G2 = np.trace(np.dot(coords1, coords1.T))
    A = np.dot(coords2.T, coords1)  # referred to as M in the original paper.
    E0 = (G1 + G2) * 0.5

    Sxx, Sxy, Sxz, Syx, Syy, Syz, Szx, Szy, Szz = A.flatten()

    Sxx2 = Sxx * Sxx
    Syy2 = Syy * Syy
    Szz2 = Szz * Szz
    Sxy2 = Sxy * Sxy
    Syz2 = Syz * Syz
    Sxz2 = Sxz * Sxz
    Syx2 = Syx * Syx
    Szy2 = Szy * Szy
    Szx2 = Szx * Szx

    SyzSzymSyySzz2 = 2.0 * (Syz * Szy - Syy * Szz)
    Sxx2Syy2Szz2Syz2Szy2 = Syy2 + Szz2 - Sxx2 + Syz2 + Szy2

    C2 = -2.0 * (Sxx2 + Syy2 + Szz2 + Sxy2 + Syx2 + Sxz2 + Szx2 + Syz2 + Szy2)
    C1 = 8.0 * (
        Sxx * Syz * Szy
        + Syy * Szx * Sxz
        + Szz * Sxy * Syx
        - Sxx * Syy * Szz
        - Syz * Szx * Sxy
        - Szy * Syx * Sxz
    )

    SxzpSzx = Sxz + Szx
    SyzpSzy = Syz + Szy
    SxypSyx = Sxy + Syx
    SyzmSzy = Syz - Szy
    SxzmSzx = Sxz - Szx
    SxymSyx = Sxy - Syx
    SxxpSyy = Sxx + Syy
    SxxmSyy = Sxx - Syy
    Sxy2Sxz2Syx2Szx2 = Sxy2 + Sxz2 - Syx2 - Szx2

    negSxzpSzx = -SxzpSzx
    negSxzmSzx = -SxzmSzx
    negSxymSyx = -SxymSyx
    SxxpSyy_p_Szz = SxxpSyy + Szz

    C0 = (
        Sxy2Sxz2Syx2Szx2 * Sxy2Sxz2Syx2Szx2
        + (Sxx2Syy2Szz2Syz2Szy2 + SyzSzymSyySzz2)
        * (Sxx2Syy2Szz2Syz2Szy2 - SyzSzymSyySzz2)
        + (negSxzpSzx * (SyzmSzy) + (SxymSyx) * (SxxmSyy - Szz))
        * (negSxzmSzx * (SyzpSzy) + (SxymSyx) * (SxxmSyy + Szz))
        + (negSxzpSzx * (SyzpSzy) - (SxypSyx) * (SxxpSyy - Szz))
        * (negSxzmSzx * (SyzmSzy) - (SxypSyx) * SxxpSyy_p_Szz)
        + (+(SxypSyx) * (SyzpSzy) + (SxzpSzx) * (SxxmSyy + Szz))
        * (negSxymSyx * (SyzmSzy) + (SxzpSzx) * SxxpSyy_p_Szz)
        + (+(SxypSyx) * (SyzmSzy) + (SxzmSzx) * (SxxmSyy - Szz))
        * (negSxymSyx * (SyzpSzy) + (SxzmSzx) * (SxxpSyy - Szz))
    )

    # Find largest root of the quaternion polynomial:
    #   f(x) = x ** 4 + c2 * x ** 2 + c1 * x + c0 = 0 (eq. 8 of Theobald et al.)
    #   f'(x) = 4 * x ** 3 + 2 * c2 * x + c1
    #
    # using Newton-Rhapson and E0 as initial guess. Liu et al. mentions 5
    # iterations are sufficient (on average) for convergence up to 1e-6
    # precision but original code writes 50, which we keep.
    nr_it = 50
    mxEigenV = E0  # starting guess (x in eqs above)
    evalprec = 1e-11  # convergence criterion
    for _ in range(nr_it):
        oldg = mxEigenV

        x2 = mxEigenV * mxEigenV
        b = (x2 + C2) * mxEigenV
        a = b + C1

        f = a * mxEigenV + C0
        f_prime = 2.0 * x2 * mxEigenV + b + a

        delta = f / (f_prime + evalprec)  # avoid division by zero
        mxEigenV = abs(mxEigenV - delta)
        if (mxEigenV - oldg) < (evalprec * mxEigenV):
            break  # convergence
    else:
        print(f"Newton-Rhapson did not converge after {nr_it} iterations")

    # The original code has a guard if minScore > 0 and rmsd < minScore, although
    # the default value of minScore is -1. For simplicity, we ignore that check.
    rmsd = (2.0 * abs(E0 - mxEigenV) / natoms) ** 0.5

    a11 = SxxpSyy + Szz - mxEigenV
    a12 = SyzmSzy
    a13 = negSxzmSzx
    a14 = SxymSyx
    a21 = SyzmSzy
    a22 = SxxmSyy - Szz - mxEigenV
    a23 = SxypSyx
    a24 = SxzpSzx
    a31 = a13
    a32 = a23
    a33 = Syy - Sxx - Szz - mxEigenV
    a34 = SyzpSzy
    a41 = a14
    a42 = a24
    a43 = a34
    a44 = Szz - SxxpSyy - mxEigenV
    a3344_4334 = a33 * a44 - a43 * a34
    a3244_4234 = a32 * a44 - a42 * a34
    a3243_4233 = a32 * a43 - a42 * a33
    a3143_4133 = a31 * a43 - a41 * a33
    a3144_4134 = a31 * a44 - a41 * a34
    a3142_4132 = a31 * a42 - a41 * a32
    q1 = a22 * a3344_4334 - a23 * a3244_4234 + a24 * a3243_4233
    q2 = -a21 * a3344_4334 + a23 * a3144_4134 - a24 * a3143_4133
    q3 = a21 * a3244_4234 - a22 * a3144_4134 + a24 * a3142_4132
    q4 = -a21 * a3243_4233 + a22 * a3143_4133 - a23 * a3142_4132

    qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4

    evecprec = 1e-6
    if qsqr < evecprec:
        q1 = a12 * a3344_4334 - a13 * a3244_4234 + a14 * a3243_4233
        q2 = -a11 * a3344_4334 + a13 * a3144_4134 - a14 * a3143_4133
        q3 = a11 * a3244_4234 - a12 * a3144_4134 + a14 * a3142_4132
        q4 = -a11 * a3243_4233 + a12 * a3143_4133 - a13 * a3142_4132
        qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4

        if qsqr < evecprec:
            a1324_1423 = a13 * a24 - a14 * a23
            a1224_1422 = a12 * a24 - a14 * a22
            a1223_1322 = a12 * a23 - a13 * a22
            a1124_1421 = a11 * a24 - a14 * a21
            a1123_1321 = a11 * a23 - a13 * a21
            a1122_1221 = a11 * a22 - a12 * a21

            q1 = a42 * a1324_1423 - a43 * a1224_1422 + a44 * a1223_1322
            q2 = -a41 * a1324_1423 + a43 * a1124_1421 - a44 * a1123_1321
            q3 = a41 * a1224_1422 - a42 * a1124_1421 + a44 * a1122_1221
            q4 = -a41 * a1223_1322 + a42 * a1123_1321 - a43 * a1122_1221
            qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4

            if qsqr < evecprec:
                q1 = a32 * a1324_1423 - a33 * a1224_1422 + a34 * a1223_1322
                q2 = -a31 * a1324_1423 + a33 * a1124_1421 - a34 * a1123_1321
                q3 = a31 * a1224_1422 - a32 * a1124_1421 + a34 * a1122_1221
                q4 = -a31 * a1223_1322 + a32 * a1123_1321 - a33 * a1122_1221
                qsqr = q1 * q1 + q2 * q2 + q3 * q3 + q4 * q4

                if qsqr < evecprec:
                    rot = np.eye(3)
                    return rmsd, rot, [q1, q2, q3, q4]

    normq = qsqr**0.5
    q1 /= normq
    q2 /= normq
    q3 /= normq
    q4 /= normq

    a2 = q1 * q1
    x2 = q2 * q2
    y2 = q3 * q3
    z2 = q4 * q4

    xy = q2 * q3
    az = q1 * q4
    zx = q4 * q2
    ay = q1 * q3
    yz = q3 * q4
    ax = q1 * q2

    rot = np.zeros((3, 3))

    rot[0][0] = a2 + x2 - y2 - z2
    rot[0][1] = 2 * (xy + az)
    rot[0][2] = 2 * (zx - ay)
    rot[1][0] = 2 * (xy - az)
    rot[1][1] = a2 - x2 + y2 - z2
    rot[1][2] = 2 * (yz + ax)
    rot[2][0] = 2 * (zx + ay)
    rot[2][1] = 2 * (yz - ax)
    rot[2][2] = a2 - x2 - y2 + z2

    return rmsd, rot, (q1, q2, q3, q4)


class QCPSuperimposer:
    """Quaternion Characteristic Polynomial (QCP) Superimposer.

    QCPSuperimposer finds the best rotation and translation to put
    two point sets on top of each other (minimizing the RMSD). This is
    eg. useful to superimposing 3D structures of proteins.

    QCP stands for Quaternion Characteristic Polynomial, which is used
    in the algorithm.

    Reference:

    Douglas L Theobald (2005), "Rapid calculation of RMSDs using a
    quaternion-based characteristic polynomial.", Acta Crystallogr
    A 61(4):478-480
    """

    def __init__(self):
        """Initialize the class."""
        self._reset_properties()

    # Private methods

    def _reset_properties(self):
        """Reset all relevant properties to None to avoid conflicts between runs."""
        self.reference_coords = None
        self.coords = None
        self.transformed_coords = None
        self.rot = None
        self.tran = None
        self.rms = None
        self.init_rms = None

    # Public methods
    def set_atoms(self, fixed, moving):
        """Prepare alignment between two atom lists.

        Put (translate/rotate) the atoms in fixed on the atoms in
        moving, in such a way that the RMSD is minimized.

        :param fixed: list of (fixed) atoms
        :param moving: list of (moving) atoms
        :type fixed,moving: [L{Atom}, L{Atom},...]
        """
        assert len(fixed) == len(moving), "Fixed and moving atom lists differ in size"

        # Grab coordinates in double precision
        fix_coord = np.array([a.get_coord() for a in fixed], dtype=np.float64)
        mov_coord = np.array([a.get_coord() for a in moving], dtype=np.float64)

        self.set(fix_coord, mov_coord)
        self.run()

        self.rms = self.get_rms()
        self.rotran = self.get_rotran()

    def apply(self, atom_list):
        """Apply the QCP rotation matrix/translation vector to a set of atoms."""
        if self.rotran is None:
            raise PDBException("No transformation has been calculated yet")

        rot, tran = self.rotran
        for atom in atom_list:
            atom.transform(rot, tran)

    # Low(er) level functions
    def set(self, reference_coords, coords):
        """Set the coordinates to be superimposed.

        coords will be put on top of reference_coords.

        - reference_coords: an NxDIM array
        - coords: an NxDIM array

        DIM is the dimension of the points, N is the number
        of points to be superimposed.
        """
        self._reset_properties()

        # store coordinates
        self.reference_coords = reference_coords
        self.coords = coords
        self._natoms, n_dim = coords.shape

        if reference_coords.shape != coords.shape:
            raise PDBException("Coordinates must have the same dimensions.")
        if n_dim != 3:
            raise PDBException("Coordinates must be Nx3 arrays.")

    def run(self):
        """Superimpose the coordinate sets."""
        if self.coords is None or self.reference_coords is None:
            raise PDBException("No coordinates set.")

        coords = self.coords.copy()
        coords_ref = self.reference_coords.copy()

        # Center Coordinates
        com_coords = np.mean(coords, axis=0)
        com_ref = np.mean(coords_ref, axis=0)

        coords -= com_coords
        coords_ref -= com_ref

        (self.rms, self.rot, _) = qcp(coords_ref, coords, self._natoms)
        self.tran = com_ref - np.dot(com_coords, self.rot)

    # Getters
    def get_transformed(self):
        """Get the transformed coordinate set."""
        if self.coords is None or self.reference_coords is None:
            raise PDBException("No coordinates set.")

        if self.rot is None:
            raise PDBException("Nothing is superimposed yet.")

        self.transformed_coords = np.dot(self.coords, self.rot) + self.tran
        return self.transformed_coords

    def get_rotran(self):
        """Return right multiplying rotation matrix and translation vector."""
        if self.rot is None:
            raise PDBException("Nothing is superimposed yet.")
        return self.rot, self.tran

    def get_init_rms(self):
        """Return the root mean square deviation of untransformed coordinates."""
        if self.coords is None:
            raise PDBException("No coordinates set yet.")

        if self.init_rms is None:
            diff = self.coords - self.reference_coords
            self.init_rms = np.sqrt(np.sum(np.sum(diff * diff, axis=1) / self._natoms))
        return self.init_rms

    def get_rms(self):
        """Root mean square deviation of superimposed coordinates."""
        if self.rms is None:
            raise PDBException("Nothing superimposed yet.")
        return self.rms
