{
 "cells": [
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "# Implementing Common Neighbor Analysis as a custom method\n",
    "\n",
    "Researchers commonly wish to implement their own custom analysis methods for particle simulations.\n",
    "Here, we show an example of how to write Common Neighbor Analysis [(Honeycutt and Andersen, J. Phys. Chem. 91, 4950)](https://pubs.acs.org/doi/abs/10.1021/j100303a014) as a custom method using `freud` and the NetworkX package.\n",
    "\n",
    "NetworkX can be installed with `pip install networkx`.\n",
    "\n",
    "First, we generate random points and determine which points share neighbors."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 1,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import defaultdict\n",
    "\n",
    "import freud\n",
    "import numpy as np"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 2,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Use a face-centered cubic (fcc) system\n",
    "box, points = freud.data.UnitCell.fcc().generate_system(4)\n",
    "aq = freud.AABBQuery(box, points)\n",
    "nl = aq.query(points, {\"num_neighbors\": 12, \"exclude_ii\": True}).toNeighborList()"
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 3,
   "metadata": {},
   "outputs": [],
   "source": [
    "# Get all sets of common neighbors.\n",
    "common_neighbors = defaultdict(list)\n",
    "for i, p in enumerate(points):\n",
    "    for j in nl.point_indices[nl.query_point_indices == i]:\n",
    "        for k in nl.point_indices[nl.query_point_indices == j]:\n",
    "            if i != k:\n",
    "                common_neighbors[(i, k)].append(j)"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Next, we use NetworkX to build graphs of common neighbors and compute the Common Neighbor Analysis signatures."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 4,
   "metadata": {},
   "outputs": [],
   "source": [
    "from collections import Counter\n",
    "\n",
    "import networkx as nx\n",
    "\n",
    "diagrams = defaultdict(list)\n",
    "particle_counts = defaultdict(Counter)\n",
    "\n",
    "for (a, b), neighbors in common_neighbors.items():\n",
    "    # Build up the graph of connections between the\n",
    "    # common neighbors of a and b.\n",
    "    g = nx.Graph()\n",
    "    for i in neighbors:\n",
    "        for j in set(nl.point_indices[nl.query_point_indices == i]).intersection(\n",
    "            neighbors\n",
    "        ):\n",
    "            g.add_edge(i, j)\n",
    "\n",
    "    # Define the identifiers for a CNA diagram:\n",
    "    # The first integer is 1 if the particles are bonded, otherwise 2\n",
    "    # The second integer is the number of shared neighbors\n",
    "    # The third integer is the number of bonds among shared neighbors\n",
    "    # The fourth integer is an index, just to ensure uniqueness of diagrams\n",
    "    diagram_type = 2 - int(b in nl.point_indices[nl.query_point_indices == a])\n",
    "    key = (diagram_type, len(neighbors), g.number_of_edges())\n",
    "    # If we've seen any neighborhood graphs with this signature,\n",
    "    # we explicitly check if the two graphs are identical to\n",
    "    # determine whether to save this one. Otherwise, we add\n",
    "    # the new graph immediately.\n",
    "    if key in diagrams:\n",
    "        isomorphs = [nx.is_isomorphic(g, h) for h in diagrams[key]]\n",
    "        if any(isomorphs):\n",
    "            idx = isomorphs.index(True)\n",
    "        else:\n",
    "            diagrams[key].append(g)\n",
    "            idx = diagrams[key].index(g)\n",
    "    else:\n",
    "        diagrams[key].append(g)\n",
    "        idx = diagrams[key].index(g)\n",
    "    cna_signature = key + (idx,)\n",
    "    particle_counts[a].update([cna_signature])"
   ]
  },
  {
   "cell_type": "markdown",
   "metadata": {},
   "source": [
    "Looking at the counts of common neighbor signatures, we see that the first particle of the fcc structure has 12 bonds with signature $(1, 4, 2, 0)$ as we expect."
   ]
  },
  {
   "cell_type": "code",
   "execution_count": 5,
   "metadata": {},
   "outputs": [
    {
     "data": {
      "text/plain": [
       "Counter({(1, 4, 2, 0): 12,\n",
       "         (2, 4, 4, 0): 6,\n",
       "         (2, 1, 0, 0): 12,\n",
       "         (2, 2, 1, 0): 24})"
      ]
     },
     "execution_count": 5,
     "metadata": {},
     "output_type": "execute_result"
    }
   ],
   "source": [
    "particle_counts[0]"
   ]
  }
 ],
 "metadata": {
  "kernelspec": {
   "display_name": "Python 3",
   "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.7.3"
  }
 },
 "nbformat": 4,
 "nbformat_minor": 2
}