From 900b09c6be9bd472da6c1e40409be72a13dbf188 Mon Sep 17 00:00:00 2001
From: ayush <ayush.suhane@gmail.com>
Date: Sat, 14 Jul 2018 14:56:55 -0700
Subject: [PATCH 1/8] added Augment functionality, tests, modified CHANGELOG

---
 package/CHANGELOG                             |   1 +
 package/MDAnalysis/lib/_augment.pyx           | 307 ++++++++++++++++++
 package/setup.py                              |   9 +-
 testsuite/MDAnalysisTests/lib/test_augment.py |  99 ++++++
 4 files changed, 415 insertions(+), 1 deletion(-)
 create mode 100644 package/MDAnalysis/lib/_augment.pyx
 create mode 100644 testsuite/MDAnalysisTests/lib/test_augment.py

diff --git a/package/CHANGELOG b/package/CHANGELOG
index 3582194490c..401adc96c89 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -20,6 +20,7 @@ The rules for this file:
 
 Enhancements
 
+  * Added augment functionality to handle periodic boundary conditions (PR #1977) 
   * Added lib.distances.capped_distance function to quickly calculate all distances
     up to a given maximum distance (PR #1941)
   * Added a rotation coordinate transformation (PR #1937)
diff --git a/package/MDAnalysis/lib/_augment.pyx b/package/MDAnalysis/lib/_augment.pyx
new file mode 100644
index 00000000000..67166d1bf11
--- /dev/null
+++ b/package/MDAnalysis/lib/_augment.pyx
@@ -0,0 +1,307 @@
+# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; -*-
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4
+#
+# MDAnalysis --- https://www.mdanalysis.org
+# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
+# (see the file AUTHORS for the full list of names)
+#
+# Released under the GNU Public Licence, v2 or any higher version
+#
+# Please cite your use of MDAnalysis in published work:
+#
+# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
+# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
+# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
+# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
+# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
+#
+# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
+# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
+# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
+#
+#
+
+import cython
+import numpy
+cimport numpy
+
+from libcpp.vector cimport vector
+from libc.math cimport sqrt
+
+__all__ = ['augment', 'undo_augment']
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def augment(float[:, ::1] coordinates, float[:, ::1] dm, float r):
+    """Calculate augmented coordinate set
+
+    Parameters
+    ----------
+    coordinates : array
+      Input coordinate array to generate duplicate images
+    dm : array
+      Real space box vectors in a matrix with shape (3, 3)
+      Vectors `[[a], [b], [c]]` are box vectors
+    r : float
+      thickness of cutoff region for duplicate image generation
+
+    Returns
+    -------
+    output : array
+      coordinates of duplicates generated due to periodic boundary conditions
+    indices : array
+      original indices of the augmented coordinates
+    """
+    cdef bint lo_x, hi_x, lo_y, hi_y, lo_z, hi_z
+    cdef int i, j, p, N
+    cdef float norm
+    cdef float shiftX[3]
+    cdef float shiftY[3]
+    cdef float shiftZ[3]
+    cdef float coord[3]
+    cdef float end[3]
+    cdef float other[3]
+    cdef float[:, ::1] reciprocal = numpy.zeros((3, 3), dtype=numpy.float32)
+    for i in range(3):
+        shiftX[i] = dm[0, i]
+        shiftY[i] = dm[1, i]
+        shiftZ[i] = dm[2, i]
+        end[i] = dm[0, i] + dm[1, i] + dm[2, i]
+    # Calculate reciprocal vectors
+    _cross(&dm[1, 0], &dm[2, 0], &reciprocal[0, 0])
+    _cross(&dm[2, 0], &dm[0, 0], &reciprocal[1, 0])
+    _cross(&dm[0, 0], &dm[1, 0], &reciprocal[2, 0])
+    # Normalize
+    for i in range(3):
+        norm = _norm(&reciprocal[i, 0])
+        for j in range(3):
+            reciprocal[i, j] = reciprocal[i, j]/norm
+
+    N = coordinates.shape[0]
+
+    cdef vector[float] output
+    cdef vector[int] indices
+
+    for i in range(0, N):
+        for j in range(3):
+            coord[j] = coordinates[i, j]
+            other[j] = end[j] - coordinates[i, j]
+        # identify the condition
+        lo_x = _dot(&coord[0], &reciprocal[0, 0]) <= r
+        hi_x = _dot(&other[0], &reciprocal[0, 0]) <= r
+        lo_y = _dot(&coord[0], &reciprocal[1, 0]) <= r
+        hi_y = _dot(&other[0], &reciprocal[1, 0]) <= r
+        lo_z = _dot(&coord[0], &reciprocal[2, 0]) <= r
+        hi_z = _dot(&other[0], &reciprocal[2, 0]) <= r
+
+        if lo_x:
+            # if X, face piece
+            for j in range(3):
+                # add to output
+                output.push_back(coord[j] + shiftX[j])
+            # keep record of which index this augmented
+            # position was created from
+            indices.push_back(i)
+
+            if lo_y:
+                # if X&Y, edge piece
+                for j in range(3):
+                    output.push_back(coord[j] + shiftX[j] + shiftY[j])
+                indices.push_back(i)
+
+                if lo_z:
+                    # if X&Y&Z, corner piece
+                    for j in range(3):
+                        output.push_back(coord[j] + shiftX[j] + shiftY[j] + shiftZ[j])
+                    indices.push_back(i)
+
+                elif hi_z:
+                    for j in range(3):
+                        output.push_back(coord[j] + shiftX[j] + shiftY[j] - shiftZ[j])
+                    indices.push_back(i)
+
+            elif hi_y:
+                for j in range(3):
+                    output.push_back(coord[j] + shiftX[j] - shiftY[j])
+                indices.push_back(i)
+
+                if lo_z:
+                    for j in range(3):
+                        output.push_back(coord[j] + shiftX[j] - shiftY[j] + shiftZ[j])
+                    indices.push_back(i)
+
+                elif hi_z:
+                    for j in range(3):
+                        output.push_back(coord[j] + shiftX[j] - shiftY[j] - shiftZ[j])
+                    indices.push_back(i)
+
+            if lo_z:
+                for j in range(3):
+                    output.push_back(coord[j] + shiftX[j] + shiftZ[j])
+                indices.push_back(i)
+
+            elif hi_z:
+                for j in range(3):
+                    output.push_back(coord[j] + shiftX[j] - shiftZ[j])
+                indices.push_back(i)
+
+        elif hi_x:
+            for j in range(3):
+                output.push_back(coord[j] - shiftX[j])
+            indices.push_back(i)
+
+            if lo_y:
+                for j in range(3):
+                    output.push_back(coord[j] - shiftX[j] + shiftY[j])
+                indices.push_back(i)
+
+                if lo_z:
+                    for j in range(3):
+                        output.push_back(coord[j] - shiftX[j] + shiftY[j] + shiftZ[j])
+                    indices.push_back(i)
+
+                elif hi_z:
+                    for j in range(3):
+                        output.push_back(coord[j] - shiftX[j] + shiftY[j] - shiftZ[j])
+                    indices.push_back(i)
+
+            elif hi_y:
+                for j in range(3):
+                    output.push_back(coord[j] - shiftX[j] - shiftY[j])
+                indices.push_back(i)
+
+                if lo_z:
+                    for j in range(3):
+                        output.push_back(coord[j] - shiftX[j] - shiftY[j] + shiftZ[j])
+                    indices.push_back(i)
+
+                elif hi_z:
+                    for j in range(3):
+                        output.push_back(coord[j] - shiftX[j] - shiftY[j] - shiftZ[j])
+                    indices.push_back(i)
+
+            if lo_z:
+                for j in range(3):
+                    output.push_back(coord[j] - shiftX[j] + shiftZ[j])
+                indices.push_back(i)
+
+            elif hi_z:
+                for j in range(3):
+                    output.push_back(coord[j] - shiftX[j] - shiftZ[j])
+                indices.push_back(i)
+
+        if lo_y:
+            for j in range(3):
+                output.push_back(coord[j] + shiftY[j])
+            indices.push_back(i)
+
+            if lo_z:
+                for j in range(3):
+                    output.push_back(coord[j] + shiftY[j] + shiftZ[j])
+                indices.push_back(i)
+
+            elif hi_z:
+                for j in range(3):
+                    output.push_back(coord[j] + shiftY[j] - shiftZ[j])
+                indices.push_back(i)
+
+        elif hi_y:
+            for j in range(3):
+                output.push_back(coord[j] - shiftY[j])
+            indices.push_back(i)
+
+            if lo_z:
+                for j in range(3):
+                    output.push_back(coord[j] - shiftY[j] + shiftZ[j])
+                indices.push_back(i)
+
+            elif hi_z:
+                for j in range(3):
+                    output.push_back(coord[j] - shiftY[j] - shiftZ[j])
+                indices.push_back(i)
+
+        if lo_z:
+            for j in range(3):
+                output.push_back(coord[j] + shiftZ[j])
+            indices.push_back(i)
+
+        elif hi_z:
+            for j in range(3):
+                output.push_back(coord[j] - shiftZ[j])
+            indices.push_back(i)
+    n = indices.size()
+    return numpy.asarray(output, dtype=numpy.float32).reshape(n, 3), numpy.asarray(indices, dtype=numpy.int32)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef float _dot(float * a, float * b):
+    """Return dot product of two sequences in range."""
+    cdef ssize_t n
+    cdef float sum1
+
+    sum1 = 0.0
+    for n in range(3):
+        sum1 += a[n] * b[n]
+    return sum1
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef void _cross(float * a, float * b, float * result):
+    """
+    Calculates the cross product between vectors
+    given by pointers a and b
+
+    Note
+    ----
+    Modifies the result array
+    """
+
+    result[0] = a[1]*b[2] - a[2]*b[1]
+    result[1] = - a[0]*b[2] + a[2]*b[0]
+    result[2] = a[0]*b[1] - a[1]*b[0]
+
+cdef float _norm(float * a):
+    """
+    Calculates the magnitude of the vector
+    """
+    cdef float result
+    cdef ssize_t n
+    result = 0.0
+    for n in range(3):
+        result += a[n]*a[n]
+    return sqrt(result)
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+def undo_augment(int[:] results, int[:] translation, int nreal):
+    """Translate augmented indices back to originals
+
+    Parameters
+    ----------
+    results : ndarray of ints
+      indices of coordinates, including "augmented" indices
+    translation : ndarray of ints
+      original indices of augmented coordinates
+    nreal : int
+      number of real coordinates, ie values in results equal or larger
+      than this need to be translated to their real counterpart
+
+    Returns
+    -------
+    results : ndarray of ints
+
+    Note
+    ----
+    Modifies the results array
+
+    """
+    cdef int N
+    N = results.shape[0]
+
+    for i in range(N):
+        if results[i] >= nreal:
+            results[i] = translation[results[i] - nreal]
+    return results
diff --git a/package/setup.py b/package/setup.py
index 4712e674d89..be8bd705b13 100755
--- a/package/setup.py
+++ b/package/setup.py
@@ -353,6 +353,13 @@ def extensions(config):
                          include_dirs=include_dirs + ['MDAnalysis/lib/include'],
                          define_macros=define_macros,
                          extra_compile_args=cpp_extra_compile_args)
+    aug = MDAExtension('MDAnalysis.lib._augment',
+                         sources=['MDAnalysis/lib/_augment' + source_suffix],
+                         language='c++',
+                         libraries=mathlib,
+                         include_dirs=include_dirs,
+                         define_macros=define_macros,
+                         extra_compile_args=cpp_extra_compile_args)
 
 
     encore_utils = MDAExtension('MDAnalysis.analysis.encore.cutils',
@@ -376,7 +383,7 @@ def extensions(config):
                               extra_compile_args=extra_compile_args)
     pre_exts = [libdcd, distances, distances_omp, qcprot,
                 transformation, libmdaxdr, util, encore_utils,
-                ap_clustering, spe_dimred, cutil]
+                ap_clustering, spe_dimred, cutil, aug]
 
     cython_generated = []
     if use_cython:
diff --git a/testsuite/MDAnalysisTests/lib/test_augment.py b/testsuite/MDAnalysisTests/lib/test_augment.py
new file mode 100644
index 00000000000..4175f5ee9eb
--- /dev/null
+++ b/testsuite/MDAnalysisTests/lib/test_augment.py
@@ -0,0 +1,99 @@
+# -*- Mode: python; tab-width: 4; indent-tabs-mode:nil; coding:utf-8 -*-
+# vim: tabstop=4 expandtab shiftwidth=4 softtabstop=4 fileencoding=utf-8
+#
+# MDAnalysis --- https://www.mdanalysis.org
+# Copyright (c) 2006-2017 The MDAnalysis Development Team and contributors
+# (see the file AUTHORS for the full list of names)
+#
+# Released under the GNU Public Licence, v2 or any higher version
+#
+# Please cite your use of MDAnalysis in published work:
+#
+# R. J. Gowers, M. Linke, J. Barnoud, T. J. E. Reddy, M. N. Melo, S. L. Seyler,
+# D. L. Dotson, J. Domanski, S. Buchoux, I. M. Kenney, and O. Beckstein.
+# MDAnalysis: A Python package for the rapid analysis of molecular dynamics
+# simulations. In S. Benthall and S. Rostrup editors, Proceedings of the 15th
+# Python in Science Conference, pages 102-109, Austin, TX, 2016. SciPy.
+#
+# N. Michaud-Agrawal, E. J. Denning, T. B. Woolf, and O. Beckstein.
+# MDAnalysis: A Toolkit for the Analysis of Molecular Dynamics Simulations.
+# J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
+
+from __future__ import absolute_import
+
+import pytest
+import numpy as np
+from numpy.testing import assert_almost_equal, assert_equal
+from itertools import product
+
+from MDAnalysis.lib._augment import augment, undo_augment
+from MDAnalysis.lib.distances import apply_PBC, transform_StoR
+from MDAnalysis.lib.mdamath import triclinic_vectors
+
+boxes = ([10, 10, 10, 90, 90, 90],  # ortho
+         [10, 10, 10, 45, 60, 90])  # tri_box
+
+# Find images for several query points,
+# here in fractional coordinates using augment
+queries = ([0.1, 0.5, 0.5],  # box face
+           [0.5, 0.5, 0.5],  # box center
+           [0.5, -0.1, 0.5],  # box face
+           [0.1, 0.1, 0.5],  # box edge
+           [0.5, -0.1, 1.1],  # box edge
+           [0.1, 0.1, 0.1],  # box vertex
+           [0.1, -0.1, 1.1],  # box vertex
+           [2.1, -3.1, 0.1]  # box vertex
+           )
+
+# Images for the previous query vectors, here in fractional coordinates.
+images = (([1.1, 0.5, 0.5],),
+          (),
+          ([0.5, -0.1, 0.5],),
+          ([1.1, 0.1, 0.5], [0.1, 1.1, 0.5], [1.1, 1.1, 0.5]),
+          ([0.5, -0.1, 0.1], [0.5, 0.9, 1.1], [0.5, -0.1, 1.1]),
+          ([1.1, 0.1, 0.1], [0.1, 1.1, 0.1], [0.1, 0.1, 1.1],
+           [0.1, 1.1, 1.1], [1.1, 1.1, 0.1], [1.1, 0.1, 1.1],
+           [1.1, 1.1, 1.1]),
+          ([1.1, 0.9, 0.1], [0.1, -0.1, 0.1], [0.1, 0.9, 1.1],
+           [0.1, -0.1, 1.1], [1.1, -0.1, 0.1], [1.1, 0.9, 1.1],
+           [1.1, -0.1, 1.1]),
+          ([1.1, 0.9, 0.1], [0.1, -0.1, 0.1], [0.1, 0.9, 1.1],
+           [0.1, -0.1, 1.1], [1.1, -0.1, 0.1], [1.1, 0.9, 1.1],
+           [1.1, -0.1, 1.1]))
+
+radius = 1.5
+
+
+@pytest.mark.parametrize('b, qres', product(boxes, zip(queries, images)))
+def test_augment(b, qres):
+    b = np.array(b, dtype=np.float32)
+    q = transform_StoR(np.array(qres[0], dtype=np.float32), b)
+    if q.shape == (3, ):
+        q = q.reshape((1, 3))
+    q = apply_PBC(q, b)
+    dm = triclinic_vectors(b)
+    aug, mapping = augment(q, dm, radius)
+    if aug.size > 0:
+        aug = np.sort(aug, axis=0)
+    else:
+        aug = list()
+    cs = transform_StoR(np.array(qres[1], dtype=np.float32), b)
+    if cs.size > 0:
+        cs = np.sort(cs, axis=0)
+    else:
+        cs = list()
+    assert_almost_equal(aug, cs, decimal=5)
+
+
+@pytest.mark.parametrize('b, qres', product(boxes, zip(queries, images)))
+def test_undoaugment(b, qres):
+    b = np.array(b, dtype=np.float32)
+    q = transform_StoR(np.array(qres[0], dtype=np.float32), b)
+    if q.shape == (3, ):
+        q = q.reshape((1, 3))
+    q = apply_PBC(q, b)
+    dm = triclinic_vectors(b)
+    aug, mapping = augment(q, dm, radius)
+    for idx, val in enumerate(aug):
+        imageid = np.asarray([len(q) + idx], dtype=np.int32)
+        assert_equal(mapping[idx], undo_augment(imageid, mapping, len(q))[0])

From 8c8edc73d2a68d9b507fa3449a3026c7a0ffdb28 Mon Sep 17 00:00:00 2001
From: ayush <ayush.suhane@gmail.com>
Date: Sat, 14 Jul 2018 15:12:20 -0700
Subject: [PATCH 2/8] Added a testcase for multiple queries

---
 testsuite/MDAnalysisTests/lib/test_augment.py | 7 +++++--
 1 file changed, 5 insertions(+), 2 deletions(-)

diff --git a/testsuite/MDAnalysisTests/lib/test_augment.py b/testsuite/MDAnalysisTests/lib/test_augment.py
index 4175f5ee9eb..687285510cd 100644
--- a/testsuite/MDAnalysisTests/lib/test_augment.py
+++ b/testsuite/MDAnalysisTests/lib/test_augment.py
@@ -42,7 +42,9 @@
            [0.5, -0.1, 1.1],  # box edge
            [0.1, 0.1, 0.1],  # box vertex
            [0.1, -0.1, 1.1],  # box vertex
-           [2.1, -3.1, 0.1]  # box vertex
+           [2.1, -3.1, 0.1],  # box vertex
+           [[0.1, 0.5, 0.5],
+            [0.5, -0.1, 0.5]] #Multiple vectors
            )
 
 # Images for the previous query vectors, here in fractional coordinates.
@@ -59,7 +61,8 @@
            [1.1, -0.1, 1.1]),
           ([1.1, 0.9, 0.1], [0.1, -0.1, 0.1], [0.1, 0.9, 1.1],
            [0.1, -0.1, 1.1], [1.1, -0.1, 0.1], [1.1, 0.9, 1.1],
-           [1.1, -0.1, 1.1]))
+           [1.1, -0.1, 1.1]),
+          ([1.1, 0.5, 0.5], [0.5, -0.1, 0.5]))
 
 radius = 1.5
 

From 0eac39b8d2a7a485965f841f9151ea7ae781edbb Mon Sep 17 00:00:00 2001
From: ayush <ayush.suhane@gmail.com>
Date: Sat, 14 Jul 2018 18:22:48 -0700
Subject: [PATCH 3/8] Modified documentation, changed the name of function

---
 package/CHANGELOG                             |  4 +-
 package/MDAnalysis/lib/_augment.pyx           | 88 ++++++++++++++-----
 testsuite/MDAnalysisTests/lib/test_augment.py | 12 ++-
 3 files changed, 75 insertions(+), 29 deletions(-)

diff --git a/package/CHANGELOG b/package/CHANGELOG
index 401adc96c89..e7b96464c60 100644
--- a/package/CHANGELOG
+++ b/package/CHANGELOG
@@ -20,7 +20,9 @@ The rules for this file:
 
 Enhancements
 
-  * Added augment functionality to handle periodic boundary conditions (PR #1977) 
+  * Added augment functionality to create relevant images of particles 
+    in the vicinity of central cell to handle periodic boundary
+    conditions (PR #1977) 
   * Added lib.distances.capped_distance function to quickly calculate all distances
     up to a given maximum distance (PR #1941)
   * Added a rotation coordinate transformation (PR #1937)
diff --git a/package/MDAnalysis/lib/_augment.pyx b/package/MDAnalysis/lib/_augment.pyx
index 67166d1bf11..693a534b2bc 100644
--- a/package/MDAnalysis/lib/_augment.pyx
+++ b/package/MDAnalysis/lib/_augment.pyx
@@ -22,38 +22,74 @@
 #
 
 import cython
-import numpy
-cimport numpy
+import numpy as np
+from .mdamath import triclinic_vectors
 
 from libcpp.vector cimport vector
 from libc.math cimport sqrt
 
-__all__ = ['augment', 'undo_augment']
+__all__ = ['augment_coordinates', 'undo_augment']
+
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def augment(float[:, ::1] coordinates, float[:, ::1] dm, float r):
-    """Calculate augmented coordinate set
+def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
+    """Calculates the relevant images of particles which are within a
+    distance 'r' from the box walls
+
+    Only the relevant images are generated i.e. for every set of
+    coordinates close to the box boundary,
+    corresponding images which are close to the opposite face and
+    outside the central cell are generated. If the particle
+    is close to more than one box walls,
+    images along the diagonals are also generated ::
+
+
+                       |  x               x
+    +---------------+  |    +---------------+
+    |               |  |    |               |
+    |               |  |    |               |
+    |               |  |    |               |
+    |             o |  |  x |             o |
+    +---------------+  |    +---------------+
+                       |
 
     Parameters
     ----------
     coordinates : array
       Input coordinate array to generate duplicate images
-    dm : array
-      Real space box vectors in a matrix with shape (3, 3)
-      Vectors `[[a], [b], [c]]` are box vectors
+      in the vicinity of the central cell. All the coordinates
+      must be within the primary unit cell.
+    box : array
+      Box dimension of shape (6, ). The dimensions must be
+      provided in the same format as returned
+      by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:
+      ``[lx, ly, lz, alpha, beta, gamma]``
     r : float
       thickness of cutoff region for duplicate image generation
 
     Returns
     -------
     output : array
-      coordinates of duplicates generated due to periodic boundary conditions
+      coordinates of duplicate(augmented) particles
     indices : array
       original indices of the augmented coordinates
+      A map which translates the indices of augmented particles
+      to their original particle index such that
+      ``indices[augmentedindex] = originalindex``
+
+    Note
+    ----
+    Output doesnot return coordinates from the initial array.
+    Use `np.concatenate(coordinates, output)` to merge particle
+    coordinates as well as their images.
+
+    .. SeeAlso:: :func:'MDAnalysis.lib._augment.undo_augment'
+
+    .. versionadded:: 0.19.0
     """
     cdef bint lo_x, hi_x, lo_y, hi_y, lo_z, hi_z
-    cdef int i, j, p, N
+    cdef int i, j, N
     cdef float norm
     cdef float shiftX[3]
     cdef float shiftY[3]
@@ -61,7 +97,11 @@ def augment(float[:, ::1] coordinates, float[:, ::1] dm, float r):
     cdef float coord[3]
     cdef float end[3]
     cdef float other[3]
-    cdef float[:, ::1] reciprocal = numpy.zeros((3, 3), dtype=numpy.float32)
+    cdef float[:, ::1] dm = np.zeros((3, 3), dtype=np.float32)
+    cdef float[:, ::1] reciprocal = np.zeros((3, 3), dtype=np.float32)
+
+    dm = triclinic_vectors(box)
+
     for i in range(3):
         shiftX[i] = dm[0, i]
         shiftY[i] = dm[1, i]
@@ -82,7 +122,7 @@ def augment(float[:, ::1] coordinates, float[:, ::1] dm, float r):
     cdef vector[float] output
     cdef vector[int] indices
 
-    for i in range(0, N):
+    for i in range(N):
         for j in range(3):
             coord[j] = coordinates[i, j]
             other[j] = end[j] - coordinates[i, j]
@@ -230,13 +270,13 @@ def augment(float[:, ::1] coordinates, float[:, ::1] dm, float r):
                 output.push_back(coord[j] - shiftZ[j])
             indices.push_back(i)
     n = indices.size()
-    return numpy.asarray(output, dtype=numpy.float32).reshape(n, 3), numpy.asarray(indices, dtype=numpy.int32)
+    return np.asarray(output, dtype=np.float32).reshape(n, 3), np.asarray(indices, dtype=np.int32)
 
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
 cdef float _dot(float * a, float * b):
-    """Return dot product of two sequences in range."""
+    """Return dot product of two 3d vectors"""
     cdef ssize_t n
     cdef float sum1
 
@@ -250,8 +290,7 @@ cdef float _dot(float * a, float * b):
 @cython.wraparound(False)
 cdef void _cross(float * a, float * b, float * result):
     """
-    Calculates the cross product between vectors
-    given by pointers a and b
+    Calculates the cross product between 3d vectors
 
     Note
     ----
@@ -277,31 +316,38 @@ cdef float _norm(float * a):
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def undo_augment(int[:] results, int[:] translation, int nreal):
-    """Translate augmented indices back to originals
+    """Translate augmented indices back to original indices
 
     Parameters
     ----------
     results : ndarray of ints
       indices of coordinates, including "augmented" indices
     translation : ndarray of ints
-      original indices of augmented coordinates
+      Map to link the augmented indices to the original particle indices
+      such that ``translation[augmentedindex] = originalindex``
     nreal : int
-      number of real coordinates, ie values in results equal or larger
+      number of real coordinates, i.e. values in results equal or larger
       than this need to be translated to their real counterpart
 
     Returns
     -------
     results : ndarray of ints
+      modified input results with all the augmented indices
+      translated to their corresponding initial original indices
 
     Note
     ----
-    Modifies the results array
+    Modifies the results array in place
+
+    .. SeeAlso:: :func:'MDAnalysis.lib._augment.augment_coordinates'
 
+    .. versionadded:: 0.19.0
     """
     cdef int N
+    cdef ssize_t i
     N = results.shape[0]
 
     for i in range(N):
         if results[i] >= nreal:
             results[i] = translation[results[i] - nreal]
-    return results
+    return np.asarray(results, dtype=np.int)
diff --git a/testsuite/MDAnalysisTests/lib/test_augment.py b/testsuite/MDAnalysisTests/lib/test_augment.py
index 687285510cd..881e275dcd1 100644
--- a/testsuite/MDAnalysisTests/lib/test_augment.py
+++ b/testsuite/MDAnalysisTests/lib/test_augment.py
@@ -26,9 +26,8 @@
 from numpy.testing import assert_almost_equal, assert_equal
 from itertools import product
 
-from MDAnalysis.lib._augment import augment, undo_augment
+from MDAnalysis.lib._augment import augment_coordinates, undo_augment
 from MDAnalysis.lib.distances import apply_PBC, transform_StoR
-from MDAnalysis.lib.mdamath import triclinic_vectors
 
 boxes = ([10, 10, 10, 90, 90, 90],  # ortho
          [10, 10, 10, 45, 60, 90])  # tri_box
@@ -44,9 +43,10 @@
            [0.1, -0.1, 1.1],  # box vertex
            [2.1, -3.1, 0.1],  # box vertex
            [[0.1, 0.5, 0.5],
-            [0.5, -0.1, 0.5]] #Multiple vectors
+            [0.5, -0.1, 0.5]]  # Multiple vectors
            )
 
+
 # Images for the previous query vectors, here in fractional coordinates.
 images = (([1.1, 0.5, 0.5],),
           (),
@@ -74,8 +74,7 @@ def test_augment(b, qres):
     if q.shape == (3, ):
         q = q.reshape((1, 3))
     q = apply_PBC(q, b)
-    dm = triclinic_vectors(b)
-    aug, mapping = augment(q, dm, radius)
+    aug, mapping = augment_coordinates(q, b, radius)
     if aug.size > 0:
         aug = np.sort(aug, axis=0)
     else:
@@ -95,8 +94,7 @@ def test_undoaugment(b, qres):
     if q.shape == (3, ):
         q = q.reshape((1, 3))
     q = apply_PBC(q, b)
-    dm = triclinic_vectors(b)
-    aug, mapping = augment(q, dm, radius)
+    aug, mapping = augment_coordinates(q, b, radius)
     for idx, val in enumerate(aug):
         imageid = np.asarray([len(q) + idx], dtype=np.int32)
         assert_equal(mapping[idx], undo_augment(imageid, mapping, len(q))[0])

From d7dd262b42abb3cb8cbe0e90dc66f681ac19a2d6 Mon Sep 17 00:00:00 2001
From: ayush <ayush.suhane@gmail.com>
Date: Sat, 14 Jul 2018 20:02:12 -0700
Subject: [PATCH 4/8] added six.moves import

---
 testsuite/MDAnalysisTests/lib/test_augment.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/testsuite/MDAnalysisTests/lib/test_augment.py b/testsuite/MDAnalysisTests/lib/test_augment.py
index 881e275dcd1..8092f133798 100644
--- a/testsuite/MDAnalysisTests/lib/test_augment.py
+++ b/testsuite/MDAnalysisTests/lib/test_augment.py
@@ -20,7 +20,7 @@
 # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
 
 from __future__ import absolute_import
-
+from six.moves import zip
 import pytest
 import numpy as np
 from numpy.testing import assert_almost_equal, assert_equal

From 776c83aa18b1db69cf21e7e5025b24034347ced1 Mon Sep 17 00:00:00 2001
From: ayush <ayush.suhane@gmail.com>
Date: Mon, 16 Jul 2018 18:27:43 -0700
Subject: [PATCH 5/8] cdef'd float instead of numpy array, changed the format
 of tests

---
 package/MDAnalysis/lib/_augment.pyx           | 34 +++----
 package/setup.py                              |  1 -
 testsuite/MDAnalysisTests/lib/test_augment.py | 96 ++++++++++---------
 3 files changed, 69 insertions(+), 62 deletions(-)

diff --git a/package/MDAnalysis/lib/_augment.pyx b/package/MDAnalysis/lib/_augment.pyx
index 693a534b2bc..5e8cf3684a9 100644
--- a/package/MDAnalysis/lib/_augment.pyx
+++ b/package/MDAnalysis/lib/_augment.pyx
@@ -97,25 +97,25 @@ def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
     cdef float coord[3]
     cdef float end[3]
     cdef float other[3]
-    cdef float[:, ::1] dm = np.zeros((3, 3), dtype=np.float32)
-    cdef float[:, ::1] reciprocal = np.zeros((3, 3), dtype=np.float32)
+    cdef float dm[3][3]
+    cdef float reciprocal[3][3]
 
     dm = triclinic_vectors(box)
 
     for i in range(3):
-        shiftX[i] = dm[0, i]
-        shiftY[i] = dm[1, i]
-        shiftZ[i] = dm[2, i]
-        end[i] = dm[0, i] + dm[1, i] + dm[2, i]
+        shiftX[i] = dm[0][i]
+        shiftY[i] = dm[1][i]
+        shiftZ[i] = dm[2][i]
+        end[i] = dm[0][i] + dm[1][i] + dm[2][i]
     # Calculate reciprocal vectors
-    _cross(&dm[1, 0], &dm[2, 0], &reciprocal[0, 0])
-    _cross(&dm[2, 0], &dm[0, 0], &reciprocal[1, 0])
-    _cross(&dm[0, 0], &dm[1, 0], &reciprocal[2, 0])
+    _cross(&dm[1][0], &dm[2][0], &reciprocal[0][0])
+    _cross(&dm[2][0], &dm[0][0], &reciprocal[1][0])
+    _cross(&dm[0][0], &dm[1][0], &reciprocal[2][0])
     # Normalize
     for i in range(3):
-        norm = _norm(&reciprocal[i, 0])
+        norm = _norm(&reciprocal[i][0])
         for j in range(3):
-            reciprocal[i, j] = reciprocal[i, j]/norm
+            reciprocal[i][j] = reciprocal[i][j]/norm
 
     N = coordinates.shape[0]
 
@@ -127,12 +127,12 @@ def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
             coord[j] = coordinates[i, j]
             other[j] = end[j] - coordinates[i, j]
         # identify the condition
-        lo_x = _dot(&coord[0], &reciprocal[0, 0]) <= r
-        hi_x = _dot(&other[0], &reciprocal[0, 0]) <= r
-        lo_y = _dot(&coord[0], &reciprocal[1, 0]) <= r
-        hi_y = _dot(&other[0], &reciprocal[1, 0]) <= r
-        lo_z = _dot(&coord[0], &reciprocal[2, 0]) <= r
-        hi_z = _dot(&other[0], &reciprocal[2, 0]) <= r
+        lo_x = _dot(&coord[0], &reciprocal[0][0]) <= r
+        hi_x = _dot(&other[0], &reciprocal[0][0]) <= r
+        lo_y = _dot(&coord[0], &reciprocal[1][0]) <= r
+        hi_y = _dot(&other[0], &reciprocal[1][0]) <= r
+        lo_z = _dot(&coord[0], &reciprocal[2][0]) <= r
+        hi_z = _dot(&other[0], &reciprocal[2][0]) <= r
 
         if lo_x:
             # if X, face piece
diff --git a/package/setup.py b/package/setup.py
index be8bd705b13..32f595e0da7 100755
--- a/package/setup.py
+++ b/package/setup.py
@@ -356,7 +356,6 @@ def extensions(config):
     aug = MDAExtension('MDAnalysis.lib._augment',
                          sources=['MDAnalysis/lib/_augment' + source_suffix],
                          language='c++',
-                         libraries=mathlib,
                          include_dirs=include_dirs,
                          define_macros=define_macros,
                          extra_compile_args=cpp_extra_compile_args)
diff --git a/testsuite/MDAnalysisTests/lib/test_augment.py b/testsuite/MDAnalysisTests/lib/test_augment.py
index 8092f133798..d15ca0404a4 100644
--- a/testsuite/MDAnalysisTests/lib/test_augment.py
+++ b/testsuite/MDAnalysisTests/lib/test_augment.py
@@ -20,57 +20,61 @@
 # J. Comput. Chem. 32 (2011), 2319--2327, doi:10.1002/jcc.21787
 
 from __future__ import absolute_import
-from six.moves import zip
 import pytest
 import numpy as np
 from numpy.testing import assert_almost_equal, assert_equal
-from itertools import product
 
 from MDAnalysis.lib._augment import augment_coordinates, undo_augment
 from MDAnalysis.lib.distances import apply_PBC, transform_StoR
 
-boxes = ([10, 10, 10, 90, 90, 90],  # ortho
-         [10, 10, 10, 45, 60, 90])  # tri_box
-
 # Find images for several query points,
-# here in fractional coordinates using augment
-queries = ([0.1, 0.5, 0.5],  # box face
-           [0.5, 0.5, 0.5],  # box center
-           [0.5, -0.1, 0.5],  # box face
-           [0.1, 0.1, 0.5],  # box edge
-           [0.5, -0.1, 1.1],  # box edge
-           [0.1, 0.1, 0.1],  # box vertex
-           [0.1, -0.1, 1.1],  # box vertex
-           [2.1, -3.1, 0.1],  # box vertex
-           [[0.1, 0.5, 0.5],
-            [0.5, -0.1, 0.5]]  # Multiple vectors
-           )
-
-
-# Images for the previous query vectors, here in fractional coordinates.
-images = (([1.1, 0.5, 0.5],),
-          (),
-          ([0.5, -0.1, 0.5],),
-          ([1.1, 0.1, 0.5], [0.1, 1.1, 0.5], [1.1, 1.1, 0.5]),
-          ([0.5, -0.1, 0.1], [0.5, 0.9, 1.1], [0.5, -0.1, 1.1]),
-          ([1.1, 0.1, 0.1], [0.1, 1.1, 0.1], [0.1, 0.1, 1.1],
-           [0.1, 1.1, 1.1], [1.1, 1.1, 0.1], [1.1, 0.1, 1.1],
-           [1.1, 1.1, 1.1]),
-          ([1.1, 0.9, 0.1], [0.1, -0.1, 0.1], [0.1, 0.9, 1.1],
-           [0.1, -0.1, 1.1], [1.1, -0.1, 0.1], [1.1, 0.9, 1.1],
-           [1.1, -0.1, 1.1]),
-          ([1.1, 0.9, 0.1], [0.1, -0.1, 0.1], [0.1, 0.9, 1.1],
-           [0.1, -0.1, 1.1], [1.1, -0.1, 0.1], [1.1, 0.9, 1.1],
-           [1.1, -0.1, 1.1]),
-          ([1.1, 0.5, 0.5], [0.5, -0.1, 0.5]))
+# here in fractional coordinates
+# Every element of qres tuple is (query, images)
+qres = (
+       ([0.1, 0.5, 0.5], [[1.1, 0.5, 0.5]]),     # box face
+       ([0.5, 0.5, 0.5], []),                    # box center
+       ([0.5, -0.1, 0.5], [[0.5, -0.1, 0.5]]),   # box face
+       ([0.1, 0.1, 0.5], [[1.1, 0.1, 0.5],
+                          [0.1, 1.1, 0.5],
+                          [1.1, 1.1, 0.5]]),     # box edge
+       ([0.5, -0.1, 1.1], [[0.5, -0.1, 0.1],
+                           [0.5, 0.9, 1.1],
+                           [0.5, -0.1, 1.1]]),   # box edge
+       ([0.1, 0.1, 0.1], [[1.1, 0.1, 0.1],
+                          [0.1, 1.1, 0.1],
+                          [0.1, 0.1, 1.1],
+                          [0.1, 1.1, 1.1],
+                          [1.1, 1.1, 0.1],
+                          [1.1, 0.1, 1.1],
+                          [1.1, 1.1, 1.1]]),     # box vertex
+       ([0.1, -0.1, 1.1], [[1.1, 0.9, 0.1],
+                           [0.1, -0.1, 0.1],
+                           [0.1, 0.9, 1.1],
+                           [0.1, -0.1, 1.1],
+                           [1.1, -0.1, 0.1],
+                           [1.1, 0.9, 1.1],
+                           [1.1, -0.1, 1.1]]),   # box vertex
+       ([2.1, -3.1, 0.1], [[1.1, 0.9, 0.1],
+                           [0.1, -0.1, 0.1],
+                           [0.1, 0.9, 1.1],
+                           [0.1, -0.1, 1.1],
+                           [1.1, -0.1, 0.1],
+                           [1.1, 0.9, 1.1],
+                           [1.1, -0.1, 1.1]]),   # box vertex
+       ([[0.1, 0.5, 0.5],
+         [0.5, -0.1, 0.5]], [[1.1, 0.5, 0.5],
+                             [0.5, -0.1, 0.5]])  # multiple queries
+       )
 
-radius = 1.5
 
-
-@pytest.mark.parametrize('b, qres', product(boxes, zip(queries, images)))
-def test_augment(b, qres):
-    b = np.array(b, dtype=np.float32)
-    q = transform_StoR(np.array(qres[0], dtype=np.float32), b)
+@pytest.mark.parametrize('b', (
+                         np.array([10, 10, 10, 90, 90, 90], dtype=np.float32),
+                         np.array([10, 10, 10, 45, 60, 90], dtype=np.float32)
+                         ))
+@pytest.mark.parametrize('q, res', qres)
+def test_augment(b, q, res):
+    radius = 1.5
+    q = transform_StoR(np.array(q, dtype=np.float32), b)
     if q.shape == (3, ):
         q = q.reshape((1, 3))
     q = apply_PBC(q, b)
@@ -79,7 +83,7 @@ def test_augment(b, qres):
         aug = np.sort(aug, axis=0)
     else:
         aug = list()
-    cs = transform_StoR(np.array(qres[1], dtype=np.float32), b)
+    cs = transform_StoR(np.array(res, dtype=np.float32), b)
     if cs.size > 0:
         cs = np.sort(cs, axis=0)
     else:
@@ -87,9 +91,13 @@ def test_augment(b, qres):
     assert_almost_equal(aug, cs, decimal=5)
 
 
-@pytest.mark.parametrize('b, qres', product(boxes, zip(queries, images)))
+@pytest.mark.parametrize('b', (
+                         np.array([10, 10, 10, 90, 90, 90], dtype=np.float32),
+                         np.array([10, 10, 10, 45, 60, 90], dtype=np.float32)
+                         ))
+@pytest.mark.parametrize('qres', qres)
 def test_undoaugment(b, qres):
-    b = np.array(b, dtype=np.float32)
+    radius = 1.5
     q = transform_StoR(np.array(qres[0], dtype=np.float32), b)
     if q.shape == (3, ):
         q = q.reshape((1, 3))

From feebb7d69887fa8d00c0c25f7fca061cb3c2702c Mon Sep 17 00:00:00 2001
From: ayush <ayush.suhane@gmail.com>
Date: Wed, 18 Jul 2018 02:57:17 -0700
Subject: [PATCH 6/8] Modified input type to undo_augment, added documentation
 stub

---
 package/MDAnalysis/lib/_augment.pyx           | 52 ++++++++++++-------
 .../documentation_pages/lib/augment.rst       |  2 +
 package/setup.py                              |  4 +-
 testsuite/MDAnalysisTests/lib/test_augment.py |  2 +-
 4 files changed, 38 insertions(+), 22 deletions(-)
 create mode 100644 package/doc/sphinx/source/documentation_pages/lib/augment.rst

diff --git a/package/MDAnalysis/lib/_augment.pyx b/package/MDAnalysis/lib/_augment.pyx
index 5e8cf3684a9..e9db62fc3e7 100644
--- a/package/MDAnalysis/lib/_augment.pyx
+++ b/package/MDAnalysis/lib/_augment.pyx
@@ -24,6 +24,7 @@
 import cython
 import numpy as np
 from .mdamath import triclinic_vectors
+cimport numpy as np
 
 from libcpp.vector cimport vector
 from libc.math cimport sqrt
@@ -37,22 +38,29 @@ def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
     """Calculates the relevant images of particles which are within a
     distance 'r' from the box walls
 
-    Only the relevant images are generated i.e. for every set of
-    coordinates close to the box boundary,
-    corresponding images which are close to the opposite face and
-    outside the central cell are generated. If the particle
-    is close to more than one box walls,
-    images along the diagonals are also generated ::
+    The algorithm works by generating explicit periodic images of
+    interior atoms. For every atom position, its distance from
+    box walls is evaluated. The distance from any respective
+    box face is calculated using a dot product of plane normal and
+    position vector of the atom. If the distance is less than a
+    specified cutoff distance, relevant periodic images are generated
+    using the box translation vectors. For instance, an atom close to
+    `XY` plane containing origin will generate a periodic image
+    outside the central cell and close to the opposite `XY` plane
+    of the box. Similarly, if the particle is close to more than
+    one box walls, images along the diagonals are also generated ::
+
+
+                           |  x               x
+        +---------------+  |    +---------------+
+        |               |  |    |               |
+        |               |  |    |               |
+        |               |  |    |               |
+        |             o |  |  x |             o |
+        +---------------+  |    +---------------+
+                           |
 
 
-                       |  x               x
-    +---------------+  |    +---------------+
-    |               |  |    |               |
-    |               |  |    |               |
-    |               |  |    |               |
-    |             o |  |  x |             o |
-    +---------------+  |    +---------------+
-                       |
 
     Parameters
     ----------
@@ -84,7 +92,10 @@ def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
     Use `np.concatenate(coordinates, output)` to merge particle
     coordinates as well as their images.
 
-    .. SeeAlso:: :func:'MDAnalysis.lib._augment.undo_augment'
+    See Also
+    --------
+    MDAnalysis.lib._augment.undo_augment
+
 
     .. versionadded:: 0.19.0
     """
@@ -270,7 +281,7 @@ def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
                 output.push_back(coord[j] - shiftZ[j])
             indices.push_back(i)
     n = indices.size()
-    return np.asarray(output, dtype=np.float32).reshape(n, 3), np.asarray(indices, dtype=np.int32)
+    return np.asarray(output, dtype=np.float32).reshape(n, 3), np.asarray(indices, dtype=np.int64)
 
 
 @cython.boundscheck(False)
@@ -315,7 +326,7 @@ cdef float _norm(float * a):
 
 @cython.boundscheck(False)
 @cython.wraparound(False)
-def undo_augment(int[:] results, int[:] translation, int nreal):
+def undo_augment(np.int64_t[:] results, np.int64_t[:] translation, int nreal):
     """Translate augmented indices back to original indices
 
     Parameters
@@ -339,7 +350,10 @@ def undo_augment(int[:] results, int[:] translation, int nreal):
     ----
     Modifies the results array in place
 
-    .. SeeAlso:: :func:'MDAnalysis.lib._augment.augment_coordinates'
+    See Also
+    --------
+    'MDAnalysis.lib._augment.augment_coordinates'
+
 
     .. versionadded:: 0.19.0
     """
@@ -350,4 +364,4 @@ def undo_augment(int[:] results, int[:] translation, int nreal):
     for i in range(N):
         if results[i] >= nreal:
             results[i] = translation[results[i] - nreal]
-    return np.asarray(results, dtype=np.int)
+    return np.asarray(results, dtype=np.int64)
diff --git a/package/doc/sphinx/source/documentation_pages/lib/augment.rst b/package/doc/sphinx/source/documentation_pages/lib/augment.rst
new file mode 100644
index 00000000000..834843d6dcb
--- /dev/null
+++ b/package/doc/sphinx/source/documentation_pages/lib/augment.rst
@@ -0,0 +1,2 @@
+.. automodule:: MDAnalysis.lib.augment
+   :members:
\ No newline at end of file
diff --git a/package/setup.py b/package/setup.py
index 32f595e0da7..6e9a7b381c0 100755
--- a/package/setup.py
+++ b/package/setup.py
@@ -353,7 +353,7 @@ def extensions(config):
                          include_dirs=include_dirs + ['MDAnalysis/lib/include'],
                          define_macros=define_macros,
                          extra_compile_args=cpp_extra_compile_args)
-    aug = MDAExtension('MDAnalysis.lib._augment',
+    augment = MDAExtension('MDAnalysis.lib._augment',
                          sources=['MDAnalysis/lib/_augment' + source_suffix],
                          language='c++',
                          include_dirs=include_dirs,
@@ -382,7 +382,7 @@ def extensions(config):
                               extra_compile_args=extra_compile_args)
     pre_exts = [libdcd, distances, distances_omp, qcprot,
                 transformation, libmdaxdr, util, encore_utils,
-                ap_clustering, spe_dimred, cutil, aug]
+                ap_clustering, spe_dimred, cutil, augment]
 
     cython_generated = []
     if use_cython:
diff --git a/testsuite/MDAnalysisTests/lib/test_augment.py b/testsuite/MDAnalysisTests/lib/test_augment.py
index d15ca0404a4..216ccfed16d 100644
--- a/testsuite/MDAnalysisTests/lib/test_augment.py
+++ b/testsuite/MDAnalysisTests/lib/test_augment.py
@@ -104,5 +104,5 @@ def test_undoaugment(b, qres):
     q = apply_PBC(q, b)
     aug, mapping = augment_coordinates(q, b, radius)
     for idx, val in enumerate(aug):
-        imageid = np.asarray([len(q) + idx], dtype=np.int32)
+        imageid = np.asarray([len(q) + idx], dtype=np.int64)
         assert_equal(mapping[idx], undo_augment(imageid, mapping, len(q))[0])

From 56c0b035dceeb3a7ee6f59828fdf2f7555d7d40b Mon Sep 17 00:00:00 2001
From: ayush <ayush.suhane@gmail.com>
Date: Thu, 19 Jul 2018 02:17:09 -0700
Subject: [PATCH 7/8] moved _dot, _cross, _norm to _cutil.pyx, modified
 document structure for _augment

---
 package/MDAnalysis/lib/_augment.pyx           | 93 +++++++------------
 package/MDAnalysis/lib/_cutil.pxd             |  3 +
 package/MDAnalysis/lib/_cutil.pyx             | 41 ++++++++
 package/MDAnalysis/lib/distances.py           |  3 +
 .../documentation_pages/lib/distances.rst     |  3 +-
 5 files changed, 83 insertions(+), 60 deletions(-)
 create mode 100644 package/MDAnalysis/lib/_cutil.pxd

diff --git a/package/MDAnalysis/lib/_augment.pyx b/package/MDAnalysis/lib/_augment.pyx
index e9db62fc3e7..a1f90631dfe 100644
--- a/package/MDAnalysis/lib/_augment.pyx
+++ b/package/MDAnalysis/lib/_augment.pyx
@@ -25,9 +25,11 @@ import cython
 import numpy as np
 from .mdamath import triclinic_vectors
 cimport numpy as np
+cimport _cutil
+from _cutil cimport _dot ,_norm, _cross
 
 from libcpp.vector cimport vector
-from libc.math cimport sqrt
+
 
 __all__ = ['augment_coordinates', 'undo_augment']
 
@@ -35,19 +37,23 @@ __all__ = ['augment_coordinates', 'undo_augment']
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
-    """Calculates the relevant images of particles which are within a
+    r"""Calculates the relevant images of particles which are within a
     distance 'r' from the box walls
 
     The algorithm works by generating explicit periodic images of
-    interior atoms. For every atom position, its distance from
-    box walls is evaluated. The distance from any respective
-    box face is calculated using a dot product of plane normal and
-    position vector of the atom. If the distance is less than a
+    interior atoms residing close to any of the six box walls. 
+    The steps involved in generating images involves
+    evaluation of reciprocal vectors for the given box vectors 
+    followed by calculation of projection distance of atom along the 
+    reciprocal vectors. If the distance is less than a
     specified cutoff distance, relevant periodic images are generated
-    using the box translation vectors. For instance, an atom close to
-    `XY` plane containing origin will generate a periodic image
+    using box translation vectors i.e. ``l[a] + m[b] + n[c]``, where 
+    ``[l, m, n]`` are the neighbouring cell indices relative to the central cell, 
+    and ``[a, b, c]`` are the box vectors. For instance, an atom close to
+    ``XY`` plane containing origin will generate a periodic image
     outside the central cell and close to the opposite `XY` plane
-    of the box. Similarly, if the particle is close to more than
+    of the box i.e. at ``0[a] + 0[b] + 1[c]``. 
+    Similarly, if the particle is close to more than
     one box walls, images along the diagonals are also generated ::
 
 
@@ -67,21 +73,21 @@ def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
     coordinates : array
       Input coordinate array to generate duplicate images
       in the vicinity of the central cell. All the coordinates
-      must be within the primary unit cell.
+      must be within the primary unit cell. (dtype = numpy.float32)
     box : array
       Box dimension of shape (6, ). The dimensions must be
       provided in the same format as returned
       by :attr:`MDAnalysis.coordinates.base.Timestep.dimensions`:
-      ``[lx, ly, lz, alpha, beta, gamma]``
+      ``[lx, ly, lz, alpha, beta, gamma]`` (dtype = numpy.float32)
     r : float
       thickness of cutoff region for duplicate image generation
 
     Returns
     -------
     output : array
-      coordinates of duplicate(augmented) particles
+      coordinates of duplicate(augmented) particles (dtype = numpy.float32)
     indices : array
-      original indices of the augmented coordinates
+      original indices of the augmented coordinates (dtype = numpy.int64)
       A map which translates the indices of augmented particles
       to their original particle index such that
       ``indices[augmentedindex] = originalindex``
@@ -89,8 +95,14 @@ def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
     Note
     ----
     Output doesnot return coordinates from the initial array.
-    Use `np.concatenate(coordinates, output)` to merge particle
-    coordinates as well as their images.
+    To merge the particles with their respective images, following operation
+    needs to be superseded after generating the images:
+    
+    .. code-block:: python
+
+            images, mapping = augment_coordinates(coordinates, box, max_cutoff)
+            all_coords = np.concatenate([coordinates, images])
+
 
     See Also
     --------
@@ -284,46 +296,6 @@ def augment_coordinates(float[:, ::1] coordinates, float[:] box, float r):
     return np.asarray(output, dtype=np.float32).reshape(n, 3), np.asarray(indices, dtype=np.int64)
 
 
-@cython.boundscheck(False)
-@cython.wraparound(False)
-cdef float _dot(float * a, float * b):
-    """Return dot product of two 3d vectors"""
-    cdef ssize_t n
-    cdef float sum1
-
-    sum1 = 0.0
-    for n in range(3):
-        sum1 += a[n] * b[n]
-    return sum1
-
-
-@cython.boundscheck(False)
-@cython.wraparound(False)
-cdef void _cross(float * a, float * b, float * result):
-    """
-    Calculates the cross product between 3d vectors
-
-    Note
-    ----
-    Modifies the result array
-    """
-
-    result[0] = a[1]*b[2] - a[2]*b[1]
-    result[1] = - a[0]*b[2] + a[2]*b[0]
-    result[2] = a[0]*b[1] - a[1]*b[0]
-
-cdef float _norm(float * a):
-    """
-    Calculates the magnitude of the vector
-    """
-    cdef float result
-    cdef ssize_t n
-    result = 0.0
-    for n in range(3):
-        result += a[n]*a[n]
-    return sqrt(result)
-
-
 @cython.boundscheck(False)
 @cython.wraparound(False)
 def undo_augment(np.int64_t[:] results, np.int64_t[:] translation, int nreal):
@@ -331,20 +303,23 @@ def undo_augment(np.int64_t[:] results, np.int64_t[:] translation, int nreal):
 
     Parameters
     ----------
-    results : ndarray of ints
-      indices of coordinates, including "augmented" indices
-    translation : ndarray of ints
+    results : numpy.ndarray
+      indices of coordinates, including "augmented" indices (dtype = numpy.int64)
+    translation : numpy.ndarray
       Map to link the augmented indices to the original particle indices
       such that ``translation[augmentedindex] = originalindex``
+      (dtype = numpy.int64)
     nreal : int
       number of real coordinates, i.e. values in results equal or larger
       than this need to be translated to their real counterpart
 
+
     Returns
     -------
-    results : ndarray of ints
+    results : numpy.ndarray
       modified input results with all the augmented indices
       translated to their corresponding initial original indices
+      (dtype = numpy.int64)
 
     Note
     ----
diff --git a/package/MDAnalysis/lib/_cutil.pxd b/package/MDAnalysis/lib/_cutil.pxd
new file mode 100644
index 00000000000..5efc82f6d5d
--- /dev/null
+++ b/package/MDAnalysis/lib/_cutil.pxd
@@ -0,0 +1,3 @@
+cdef float _dot(float *, float *)
+cdef void _cross(float *, float *, float *)
+cdef float _norm(float *)
\ No newline at end of file
diff --git a/package/MDAnalysis/lib/_cutil.pyx b/package/MDAnalysis/lib/_cutil.pyx
index 0463cf15cef..ef1b5f00b3e 100644
--- a/package/MDAnalysis/lib/_cutil.pyx
+++ b/package/MDAnalysis/lib/_cutil.pyx
@@ -24,6 +24,7 @@
 import cython
 import numpy as np
 cimport numpy as np
+from libc.math cimport sqrt
 
 from MDAnalysis import NoDataError
 
@@ -262,3 +263,43 @@ def make_whole(atomgroup, reference_atom=None):
         raise ValueError("AtomGroup was not contiguous from bonds, process failed")
     else:
         atomgroup.positions = newpos
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef float _dot(float * a, float * b):
+    """Return dot product of two 3d vectors"""
+    cdef ssize_t n
+    cdef float sum1
+
+    sum1 = 0.0
+    for n in range(3):
+        sum1 += a[n] * b[n]
+    return sum1
+
+
+@cython.boundscheck(False)
+@cython.wraparound(False)
+cdef void _cross(float * a, float * b, float * result):
+    """
+    Calculates the cross product between 3d vectors
+
+    Note
+    ----
+    Modifies the result array
+    """
+
+    result[0] = a[1]*b[2] - a[2]*b[1]
+    result[1] = - a[0]*b[2] + a[2]*b[0]
+    result[2] = a[0]*b[1] - a[1]*b[0]
+
+cdef float _norm(float * a):
+    """
+    Calculates the magnitude of the vector
+    """
+    cdef float result
+    cdef ssize_t n
+    result = 0.0
+    for n in range(3):
+        result += a[n]*a[n]
+    return sqrt(result)
\ No newline at end of file
diff --git a/package/MDAnalysis/lib/distances.py b/package/MDAnalysis/lib/distances.py
index 790e540ed08..907bbd44cbf 100644
--- a/package/MDAnalysis/lib/distances.py
+++ b/package/MDAnalysis/lib/distances.py
@@ -61,6 +61,8 @@
 .. autofunction:: apply_PBC(coordinates, box [, backend])
 .. autofunction:: transform_RtoS(coordinates, box [, backend])
 .. autofunction:: transform_StoR(coordinates, box [,backend])
+.. autofunction:: MDAnalysis.lib._augment.augment_coordinates(coordinates, box, radius)
+.. autofunction:: MDAnalysis.lib._augment.undo_augment(indices, translation, nreal)
 
 """
 from __future__ import division, absolute_import
@@ -70,6 +72,7 @@
 from numpy.lib.utils import deprecate
 
 from .mdamath import triclinic_vectors, triclinic_box
+from ._augment import augment_coordinates, undo_augment
 
 
 
diff --git a/package/doc/sphinx/source/documentation_pages/lib/distances.rst b/package/doc/sphinx/source/documentation_pages/lib/distances.rst
index 803d8c7fd0a..8d89a94b09c 100644
--- a/package/doc/sphinx/source/documentation_pages/lib/distances.rst
+++ b/package/doc/sphinx/source/documentation_pages/lib/distances.rst
@@ -2,6 +2,7 @@
    :members:
 
 
+
 Low-level modules for :mod:`MDAnalysis.lib.distances`
 =====================================================
 
@@ -15,4 +16,4 @@ the OpenMP-enable functions.
    :maxdepth: 1
 
    c_distances
-   c_distances_openmp
+   c_distances_openmp
\ No newline at end of file

From 7ef48ef8f36ca795974d8be8f3dc9825c12c0394 Mon Sep 17 00:00:00 2001
From: ayush <ayush.suhane@gmail.com>
Date: Thu, 19 Jul 2018 02:30:06 -0700
Subject: [PATCH 8/8] removed augment.rst

---
 package/doc/sphinx/source/documentation_pages/lib/augment.rst | 2 --
 1 file changed, 2 deletions(-)
 delete mode 100644 package/doc/sphinx/source/documentation_pages/lib/augment.rst

diff --git a/package/doc/sphinx/source/documentation_pages/lib/augment.rst b/package/doc/sphinx/source/documentation_pages/lib/augment.rst
deleted file mode 100644
index 834843d6dcb..00000000000
--- a/package/doc/sphinx/source/documentation_pages/lib/augment.rst
+++ /dev/null
@@ -1,2 +0,0 @@
-.. automodule:: MDAnalysis.lib.augment
-   :members:
\ No newline at end of file