{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "### Atom-Atom Mapping\n",
    "\n",
    "Given two molecular structures, it is important to identify atoms that are chemically similar. This a commonly used in 3D-QSAR pharmacore analysis, substructure searching, metabolic pathway identification, and chemical machine learning. This problem can be formulated as a 2-sided permutation Procrustes with single transformation.\n",
    "\n",
    "> **Permutation Procrustes 2-Sided with Single-Transformation**\n",
    ">\n",
    "> Given matrix $\\mathbf{A}_{n \\times n}$ and a reference $\\mathbf{B}_{n \\times n}$, find a permutation of rows/columns of $\\mathbf{A}_{n \\times n}$ that makes it as close as possible to $\\mathbf{B}_{n \\times n}$, i.e.,\n",
    ">\n",
    "\\begin{equation}\n",
    "   \\underbrace{\\text{min}}_{\\left\\{\\mathbf{P} \\left| {p_{ij} \\in \\{0, 1\\}\n",
    "               \\atop \\sum_{i=1}^n p_{ij} = \\sum_{j=1}^n p_{ij} = 1} \\right. \\right\\}}\n",
    "               \\|\\mathbf{P}^\\dagger \\mathbf{A} \\mathbf{P} - \\mathbf{B}\\|_{F}^2 \\\\\n",
    "\\end{equation}\n",
    "\n",
    "\n",
    "In the code block below, we use the `procrustes` library to map atoms of *but-1-en-3-yne* (molecule **A**) and *3,3-dimethylpent-1-en-4-yne* (molecule **B**) in **Fig. (i)**. Based on our chemical intuition, we can tell that the triple and double bonds of the molecules \"match\"; however, simple (geometric) molecular alignment based on three-dimensional coordinates does not identify that. The key step is defining a representation that contains bonding information before applying permutation Procrustes to match atoms. \n",
    "\n",
    "- **Fig. (ii):** Inspired by graph theory, we represent each molecule with an \"adjacency\" matrix where the diagonal elements are the atomic numbers and the off-diagonal elements are the bond orders. This results in matrices $\\mathbf{A} \\in \\mathbb{R}^{4 \\times 4}$ and $\\mathbf{B} \\in \\mathbb{R}^{7 \\times 7}$. Note that the permutation Procrustes requires the two matrices to be of the same size, so the smaller matrix $\\mathbf{A}$ is padded with zero rows and columns to have same shape as matrix $\\mathbf{B}$.\n",
    "\n",
    "- **Fig. (iii):** The mapping between atoms can be also directly deduced from the optimal permutation matrix $\\mathbf{P}$. Specifically, the transformed matrix $\\mathbf{P^{\\top}AP}$ should be compared to matrix $\\mathbf{B}$ to identify the matching atoms; the zero rows/columns in $\\mathbf{A}$ (colored in blue) correspond to atoms in $\\mathbf{B}$ for which there are no corresponding atoms. "
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "![Fig. 1. Atom-atom Mapping with Two-sided Permutation Procrustes](notebook_data/atom_atom_mapping/atom_atom_mapping.png \"Fig. 1 Atom-atom Mapping with Two-sided Permutation Procrustes\")"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [
    {
     "name": "stdout",
     "output_type": "stream",
     "text": [
      "Permutation Matrix:\n",
      " [[1. 0. 0. 0. 0. 0. 0.]\n",
      " [0. 1. 0. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 1. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 1. 0. 0.]\n",
      " [0. 0. 1. 0. 0. 0. 0.]\n",
      " [0. 0. 0. 0. 0. 1. 0.]\n",
      " [0. 0. 0. 0. 0. 0. 1.]]\n",
      "\n",
      "Transformed A: \n",
      " [[6 3 0 0 0 0 0]\n",
      " [3 6 0 1 0 0 0]\n",
      " [0 0 0 0 0 0 0]\n",
      " [0 1 0 6 2 0 0]\n",
      " [0 0 0 2 6 0 0]\n",
      " [0 0 0 0 0 0 0]\n",
      " [0 0 0 0 0 0 0]]\n",
      "\n",
      "Compare to Original (padded) B:\n",
      " [[6 3 0 0 0 0 0]\n",
      " [3 6 1 0 0 0 0]\n",
      " [0 1 6 1 0 1 1]\n",
      " [0 0 1 6 2 0 0]\n",
      " [0 0 0 2 6 0 0]\n",
      " [0 0 1 0 0 6 0]\n",
      " [0 0 1 0 0 0 6]]\n",
      "\n",
      "Procrustes Error: 118.0\n"
     ]
    }
   ],
   "source": [
    "# atom-atom mapping with 2-sided permutation procrustes (with single transformation)\n",
    "\n",
    "import numpy as np\n",
    "from procrustes import permutation_2sided\n",
    "\n",
    "# Define molecule A representing but‐1‐en‐3‐yne\n",
    "A = np.array([[6, 3, 0, 0],\n",
    "              [3, 6, 1, 0],\n",
    "              [0, 1, 6, 2],\n",
    "              [0, 0, 2, 6]])\n",
    "\n",
    "# Define molecule B representing 3,3‐dimethylpent‐1‐en‐4‐yne\n",
    "B = np.array([[6, 3, 0, 0, 0, 0, 0],\n",
    "              [3, 6, 1, 0, 0, 0, 0],\n",
    "              [0, 1, 6, 1, 0, 1, 1],\n",
    "              [0, 0, 1, 6, 2, 0, 0],\n",
    "              [0, 0, 0, 2, 6, 0, 0],\n",
    "              [0, 0, 1, 0, 0, 6, 0],\n",
    "              [0, 0, 1, 0, 0, 0, 6]])\n",
    "\n",
    "# Two-sided permutation Procrustes (with single transformation)\n",
    "result = permutation_2sided(A, B, method=\"approx-normal1\", single=True, pad=True)\n",
    "\n",
    "# Compute the transformed molecule A using transformation matrix P\n",
    "P = result.t\n",
    "new_A = np.dot(P.T, np.dot(result.new_a, P)).astype(int)\n",
    "\n",
    "print(\"Permutation Matrix:\\n\", P)\n",
    "print(\"\\nTransformed A: \\n\", new_A)\n",
    "print(\"\\nCompare to Original (padded) B:\\n\", result.new_b)\n",
    "print(\"\\nProcrustes Error:\", result.error)"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3 (ipykernel)",
   "language": "python",
   "name": "python3"
  },
  "language_info": {
   "codemirror_mode": {
    "name": "ipython",
    "version": 3
   },
   "file_extension": ".py",
   "mimetype": "text/x-python",
   "name": "python",
   "nbconvert_exporter": "python",
   "pygments_lexer": "ipython3",
   "version": "3.8.12"
  },
  "toc": {
   "base_numbering": 1,
   "nav_menu": {},
   "number_sections": true,
   "sideBar": true,
   "skip_h1_title": false,
   "title_cell": "Table of Contents",
   "title_sidebar": "Contents",
   "toc_cell": false,
   "toc_position": {},
   "toc_section_display": true,
   "toc_window_display": false
  }
 },
 "nbformat": 4,
 "nbformat_minor": 4
}