added some comments and created lifted directory

pyqc/lift.py

@@ -26,6 +26,22 @@ def __init__(self,
                  lattice: np.ndarray,
                  is2D: bool
                  ):
+        """
+        Parameters
+        ----------
+        box: freud.box.Box object
+        pos: N x 3 numpy array of float64
+        query_args: dictionary containing arguments for freud.AABBQuery. Default is None
+        hd_lattice: N x M numpy array of ints, where M is the hyperlattice dimensions
+        avg_a: float64
+        lifted_pos: N_0 x 3 numpy array of float 64 (N_0 <= N)
+        a: float64
+        projected_dim: unsigned int, number of hyperlattice dimensions (M)
+        parallel_basis: M x 3 numpy array of float 64
+        perp_basis: M x 3 numpy array of float 64
+        lattice: M x M numpy array of ints
+        is2D: bool
+        """
 
         self._a = a
         self._projected_dim = projected_dim
@@ -95,24 +111,40 @@ def is2D(self):
     def compute(self, origin: np.ndarray,
                 tolerance: float = 0.05,
                 verbose: bool = False) -> None:
-
         origin = np.atleast_2d(origin)
         pos3d = np.copy(self._pos)
-
+        """
+        origin: 1 x 3 numpy array of float64 [[x, y, z]]
+        tolerance: float
+        vebose: bool
+        """
         nl = freud.locality.AABBQuery(
             self._box, pos3d).query(
             pos3d, self._query_args).toNeighborList()
+        # -- freud Neighborlist Object
+        """To get an array:
+        nl_array = nl[:]
+        # -- nl_array: K x 2 array of unsigned ints, K = len(nl_array)
+        """
 
         # Sort positions based on distance from origin to make lifting more efficient
         r0j = self._box.compute_all_distances(pos3d, origin).flatten()
+        # -- r0j: N x 1 numpy array of floats
+        # -- N = len(pos3d)
         dist_sorted = np.argsort(r0j)
+        # -- dist_sorted: N x 1 numpy array of unsigned ints
 
         pos = np.copy(pos3d)
+        # -- pos: N x 3 numpy array of floats, the size of this matrix changes
         if self.is2D:
-            pos[:, -1] = 0.
+            pos[:, -1] = 0. # confused about what this means: why -1? (last column)
+            # including ... means the last index
 
         hd_coord = np.zeros([len(pos), self._projected_dim])
+        # -- hd_coord: N x M numpy array of floats, where N = len(pos) and M = self._projected_dim
+        # ----- These should honestly be ints but somewhere in Python land, some of my int arrays turned into floats,
         mag_basis = self._box.compute_all_distances(self._parallel_basis, [[0, 0, 0]]).flatten()
+        # -- mag_basis: M x 1 numpy array of floats
 
         # For average tile length calc, store sum of tiles for each lifted case
         a_avg = 0.0
@@ -124,30 +156,90 @@ def compute(self, origin: np.ndarray,
         # Presort nl to reduce in1d calls
         # nl_red, nl_construct = np.unique(nl[:, 0], return_inverse=True)
 
+    # don't really completely understand what this section is doing
+        # getting index of particles already lifted
         ind_lift = np.empty(len(pos)).astype(int)
+        # outer shell of particles lifting from to find neighbors
+        # -- ind_lift: N x 1 array, int
         ind_ref = np.atleast_1d(dist_sorted[n_tried]).astype(int)
+        # -- ind_ref: N x 1 array, int (can use unsigned int)
         ind_lift[n_lifted] = dist_sorted[n_tried]
         # ind_red = np.arange(len(pos))
         # ind_red = ind_red[ind_red != dist_sorted[n_tried]]
         pbar = tqdm(total=len(pos))
         while n_tried < len(pos):
             pbar.update(n_lifted + n_tried)
+            # in freud: indexes of list of neighbors, calculated everything at once
             nl_mod = nl[np.in1d(nl[:, 0], ind_ref)]
+            # -- nl_mod: L x 1 array of unsigned ints
+            # --- nl[:, 0] are the indicies of reference points and nl[:, 1] are the indices of their neighbors
+            """
+            z = np.in1d(x, y) 
+            x: a 1D array
+            y: a 1D array
+            z: a 1D array of booleans with size len(x)
+            Does the same thing as the following code (but way faster, loops in Python are slow):
+
+            z = np.zeros_like(x).astype(bool)  # Makes a 1D array of False with same dimensions as x
+            for n in range(len(x)):
+                for _y in y:
+                    if x[n] == _y:
+                        z[n] = True
+            """
             nl_mod = nl_mod[np.in1d(nl_mod[:, 1], ind_lift, invert=True)]
+            # computing vectors
+            # -- nl_mod: 1D array of unsigned ints, size is J = len(nl_mod)
+            # ---- Removes particles that have already been lifted from the list of candidates
+            """
+            z = np.in1d(x, y, invert=True) 
+            x: a 1D array
+            y: a 1D array
+            z: a 1D array of booleans with size len(x)
+            Does the same thing as the following code (but way faster, loops in Python are slow):
+
+            z = np.ones_like(x).astype(bool)  # Makes a 1D array of True with same dimensions as x
+            for n in range(len(x)):
+                for _y in y:
+                    if x[n] == _y:
+                        z[n] = False
+            """
             vij = self._box.wrap(pos[nl_mod[:, 1]] - pos[nl_mod[:, 0]])
+            # -- vij = J x 3 array of float32
+            # --- need to wrap vectors because the box is periodic
+            # --- Source code https://github.com/glotzerlab/freud/blob/3b9b38778e248ae5df6e0d5fd502924e62e06fe7/cpp/box/Box.h#L311-L333
             vij /= self._box.compute_all_distances(vij, [[0, 0, 0]])
+            # -- vij = J x 3 array of float32
+            # --- https://github.com/glotzerlab/freud/blob/3b9b38778e248ae5df6e0d5fd502924e62e06fe7/cpp/box/Box.h#L456-L469
+
             # Find particles which align with most with basis
             dot_prod = vij @ self._parallel_basis.T
-            dot_max = np.amax(abs(dot_prod), axis=1)
-            basis_ind = np.argmax(abs(dot_prod), axis=1)
+            # -- dot_prod = J x M array of floats, where M = self._projected_dim
+            dot_max = np.amax(abs(dot_prod), axis=1) # Returns maximum value of dot_prod for each row
+            # -- dot_max: J x 1 array of floats
+            # amax gives you the actual value
+            basis_ind = np.argmax(abs(dot_prod), axis=1) # Returns column index of the maximum value of dot_prod, per row
+            # -- basis_ind: J x 1 array of unsigned ints
+            # argmax tells you where it is
 
             # Check that the angle of alignment is within tolerance level
             align_ind = np.where(dot_max >= (1-tolerance)*mag_basis[basis_ind])[0]
+            # -- align_ind: 1D array of unsigned unsigned ints
 
             # Get indicies of aligned positions for lifting
             if len(align_ind) > 0:
                 # Prevent double lifting (if particles 1 and 2 both lift the same particle)
                 lifted_pos_ind, nl_ind, counts = np.unique(nl_mod[align_ind, 1], return_index=True, return_counts=True)
+                # -- lifted_pos_ind: 1D array of unsigned ints
+                # -- nl_ind: 1D array of unsigned ints
+                # -- counts: 1D array of unsigned ints
+                """
+                Example:
+                nl_mod = [[0, 12], [0, 3], [1, 4], [1, 12]]
+                lifted_pos_ind, nl_ind, counts = np.unique(nl_mod[:, 1], return_index=True, return_counts=True)
+                lifted_pos_ind = [3, 4, 12]  This doesn't need to be in order, but np.unique sorts
+                nl_ind = [1, 2, 0]
+                counts = [1, 1, 2]
+                """
                 ind_repeat = lifted_pos_ind[np.where(counts > 1)[0]]
                 for ind in ind_repeat:
                     nl_repeat = align_ind[nl_mod[align_ind, 1] == ind]
@@ -190,7 +282,7 @@ def compute(self, origin: np.ndarray,
                 # plt.gca().set_aspect('equal')
                 # plt.show()
 
-                if verbose:
+                if verbose: # what does verbose stand for here? prints out where you are in the printing
                     logging.info("%i of %i particles lifted" % (n_lifted, len(pos)))
 
             elif n_lifted <= 1:

pyqc/qc.py

@@ -4,11 +4,16 @@
 from .lift import Lift
 from .utils import make_unit, basis_values, alignment
 
+# decimal expansion of golden ratio
 _tau = 0.5 * (1 + np.sqrt(5))
 
 
 class QC:
+    """
+    Specify information about the type of quasicrystal.
+    """
 
+    # constructor
     def __init__(
             self,
             a: float,
@@ -62,18 +67,35 @@ def ref_vector(self, value: np.ndarray):
 
     @property
     def projected_dim(self) -> int:
+        """
+        Projected dimension -- either 2D or 3D.
+        :return:
+        """
         return self._projected_dim
 
     @property
     def parallel_basis(self) -> np.ndarray:
+        """
+        Set of vectors that will help you travel to each point in the lattice.
+        :return:
+        """
         return self._parallel_basis
 
     @property
     def perp_basis(self) -> np.ndarray:
+        """
+
+        :return:
+        """
         return self._perp_basis
 
     @property
     def lattice(self) -> np.ndarray:
+        """
+        Higher dimensional lattice
+        Usually an identity matrix with the same dimensions as the hyperspace.
+        :return:
+        """
         return self._lattice
 
     @lattice.setter
@@ -98,6 +120,15 @@ def axial(
             ref_vector: np.ndarray or list = None,
             perp_vector: np.ndarray or list = None
     ):
+        """
+        Class of quasicrystals (quasiperiodic in 2D, periodic in 1D)
+        :param a:
+        :param theta_par:
+        :param theta_perp:
+        :param ref_vector:
+        :param perp_vector:
+        :return:
+        """
         if ref_vector is not None and perp_vector is not None:
             raise ValueError
         else:
@@ -131,6 +162,14 @@ def penrose(
             ref_vector: np.ndarray or list = None,
             perp_vector: np.ndarray or list = None
     ):
+        """
+        Penrose tiling.
+        :param a: Bond length in system
+        :param ref_vector: Reference vector for basis set. The basis set will be rotated so that
+        the first vector in the basis lines up with the reference vector.
+        :param perp_vector:
+        :return:
+        """
         theta_par = 2*np.pi/5
         theta_perp = 2*theta_par
         pen = cls.axial(a, theta_par, theta_perp,

requirements-testing.txt

@@ -1,7 +1,10 @@
 numpy~=1.20.0rc2
 rowan~=1.3.0.post1
-freud~=0.1.2
+freud-analysis~=2.5.1
 scipy~=1.6.0
 sklearn~=0.0
 scikit-learn~=0.24.1
 setuptools~=49.6.0
+tqdm~=4.56.0
+matplotlib~=3.3.3
+gsd~=2.4.0

requirements.txt

@@ -1,10 +1,9 @@
 numpy~=1.20.0rc2
 rowan~=1.3.0.post1
-freud~=0.1.2
+freud-analysis~=2.5.1
 scipy~=1.6.0
 sklearn~=0.0
 scikit-learn~=0.24.1
 setuptools~=49.6.0
 tqdm~=4.56.0
-gsd~=2.4.0
 matplotlib~=3.3.3

test/test.py

@@ -1,4 +1,11 @@
 #!/usr/bin/env python3
+
+"""
+import sys
+sys.path
+sys.path.append('../')
+"""
+
 import pyqc
 
 import numpy as np
@@ -14,16 +21,19 @@
 
 def main():
     # all_files = glob.glob('*.gsd', recursive=True)
-    all_files = glob.glob('dec_swap*.gsd', recursive=True)
+    # glob.glob: gives a list of all the files that match the string pattern specified
+    all_files = glob.glob('*.gsd', recursive=True)
     ignore_files = glob.glob('lifted/*.gsd', recursive=True)
-    ignore_files.extend(glob.glob('test*.gsd', recursive=True))
+    ignore_files.extend(glob.glob('dec_swap*.gsd', recursive=True))
+    ignore_files.append('dodec.gsd')
+    ignore_files.append('ideal_tiling.gsd')
     test_files = []
     for f in all_files:
         if f not in ignore_files:
             test_files.append(f)
-    test_files.append('ideal_tiling.gsd')
-    # test_files = ['ideal_tiling.gsd', 'dodec.gsd',
-    #               'PhiPrime.gsd', '21.gsd', '32.gsd', '53.gsd']  #kwangs - 1.22e-4
+    # test_files.append('ideal_tiling.gsd')
+    # test_files = ['test*.gsd', 'PhiPrime.gsd',
+    #               '21.gsd', '32.gsd', '53.gsd']  #kwangs - 1.22e-4
     ps = []
     ps_old = []
     for f in test_files:
@@ -43,6 +53,7 @@ def main():
             }
             qc = pyqc.QC.icosahedral(a=1.8)
             dr = 0.6
+            """
         elif 'dodec.gsd' == f:
             a = 1.1
             query_args = {
@@ -61,6 +72,7 @@ def main():
                 'exclude_ii': True
             }
             qc = pyqc.QC.decagonal(a=1.1)
+            """
         else:
             query_args = {
                 'mode': 'ball',
@@ -94,6 +106,8 @@ def main():
             pass
 
         # local = ps_.compute_local()
+        if os.path.isdir('lifted'):
+            os.mkdir('lifted')
 
         f_proj = gsd.hoomd.open(f'lifted/lifted_{f}', mode='wb')