Reverting HexOrderParameter changes to master.

cpp/order/HexOrderParameter.cc

@@ -16,9 +16,7 @@ using namespace tbb;
 
 namespace freud { namespace order {
 
-HexOrderParameter::HexOrderParameter(float rmax,
-                                     unsigned int k,
-                                     unsigned int n)
+HexOrderParameter::HexOrderParameter(float rmax, float k, unsigned int n)
     : m_box(box::Box()), m_k(k), m_Np(0)
     {
     }
@@ -27,55 +25,50 @@ HexOrderParameter::~HexOrderParameter()
     {
     }
 
-void HexOrderParameter::compute(box::Box& box,
-                                const freud::locality::NeighborList *nlist,
-                                const vec3<float> *points,
-                                unsigned int Np)
+void HexOrderParameter::compute(box::Box& box, const freud::locality::NeighborList *nlist,
+                                const vec3<float> *points, unsigned int Np)
     {
-    // Compute the cell list
+    // compute the cell list
     m_box = box;
 
     nlist->validate(Np, Np);
     const size_t *neighbor_list(nlist->getNeighbors());
 
-    // Reallocate the output array if it is not the right size
+    // reallocate the output array if it is not the right size
     if (Np != m_Np)
         {
-        m_psi_array = std::shared_ptr<complex<float> >(new complex<float> [Np],
-                std::default_delete<complex<float>[]>());
+        m_psi_array = std::shared_ptr<complex<float> >(new complex<float> [Np], std::default_delete<complex<float>[]>());
         }
 
-    // Compute the order parameter
+    // compute the order parameter
     parallel_for(blocked_range<size_t>(0,Np),
         [=] (const blocked_range<size_t>& r)
         {
         size_t bond(nlist->find_first_index(r.begin()));
-        for (size_t i=r.begin(); i!=r.end(); ++i)
+        for(size_t i=r.begin(); i!=r.end(); ++i)
             {
             m_psi_array.get()[i] = 0;
             vec3<float> ref = points[i];
 
-            for (; bond < nlist->getNumBonds() && neighbor_list[2*bond] == i;
-                ++bond)
+            for(; bond < nlist->getNumBonds() && neighbor_list[2*bond] == i; ++bond)
                 {
                 const size_t j(neighbor_list[2*bond + 1]);
 
-                // Compute r between the two particles
+                //compute r between the two particles
                 vec3<float> delta = m_box.wrap(points[j] - ref);
 
                 float rsq = dot(delta, delta);
                 if (rsq > 1e-6)
                     {
-                    // Compute psi for neighboring particle
-                    // (only constructed for 2d)
+                    //compute psi for neighboring particle(only constructed for 2d)
                     float psi_ij = atan2f(delta.y, delta.x);
-                    m_psi_array.get()[i] += exp(complex<float>(0, m_k*psi_ij));
+                    m_psi_array.get()[i] += exp(complex<float>(0,m_k*psi_ij));
                     }
                 }
             m_psi_array.get()[i] /= complex<float>(m_k);
             }
         });
-    // Save the last computed number of particles
+    // save the last computed number of particles
     m_Np = Np;
     }
 

cpp/order/HexOrderParameter.h

@@ -35,7 +35,7 @@ class HexOrderParameter
     {
     public:
         //! Constructor
-        HexOrderParameter(float rmax, unsigned int k=6, unsigned int n=0);
+        HexOrderParameter(float rmax, float k=6, unsigned int n=0);
 
         //! Destructor
         ~HexOrderParameter();
@@ -63,14 +63,14 @@ class HexOrderParameter
             return m_Np;
             }
 
-        unsigned int getK()
+        float getK()
             {
             return m_k;
             }
 
     private:
         box::Box m_box;            //!< Simulation box where the particles belong
-        unsigned int m_k;          //!< Multiplier in the exponent
+        float m_k;                 //!< Multiplier in the exponent
         unsigned int m_Np;         //!< Last number of points computed
 
         std::shared_ptr< std::complex<float> > m_psi_array;  //!< psi array computed

cpp/order/wigner3j.h

@@ -12,36 +12,8 @@
 
 using namespace std;
 
-vector<float> getWigner3j(unsigned int l);
-
 // All wigner3j coefficients created using sympy
-/*
-from sympy.physics.wigner import wigner_3j
-import scipy.io as sio
-import numpy as np
-import json
-
-wigner = {}
-for i in range(10):
-    counter = 0
-    l = (i+1)*2
-    print('Computing l={}'.format(l))
-    for u1 in range(2*l+1):
-        for u2 in range(max(l-u1,0),min(2*l+1,3*l-u1+1)):
-            counter += 1
-    W_l = np.zeros(counter, dtype=np.float64)
-
-    j = 0
-    for u1 in range(2*l+1):
-        for u2 in range(max(l-u1,0),min(2*l+1,3*l-u1+1)):
-            u3 = 3*l-u2-u1
-            W_l[j] = float(wigner_3j(l,l,l,u1-l,u2-l,u3-l))
-            j += 1
-    wigner['l'+str(l)] = W_l.tolist()
-
-with open('wigner3j.json', 'w') as jsonfile:
-    json.dump(wigner, jsonfile, sort_keys=True, indent=4)
-sio.savemat('wigner3j.mat',wigner)'''
-*/
+// http://www.sympy.org/
+vector<float> getWigner3j(unsigned int l);
 
 #endif

freud/_order.pxd

@@ -55,7 +55,7 @@ cdef extern from "CubaticOrderParameter.h" namespace "freud::order":
 
 cdef extern from "HexOrderParameter.h" namespace "freud::order":
     cdef cppclass HexOrderParameter:
-        HexOrderParameter(float, unsigned int, unsigned int)
+        HexOrderParameter(float, float, unsigned int)
         const box.Box &getBox() const
         void compute(box.Box &,
                      const freud._locality.NeighborList*,
@@ -64,7 +64,7 @@ cdef extern from "HexOrderParameter.h" namespace "freud::order":
         # unsure how to pass back the std::complex, but this seems to compile...
         shared_array[float complex] getPsi()
         unsigned int getNP()
-        unsigned int getK()
+        float getK()
 
 cdef extern from "LocalDescriptors.h" namespace "freud::order":
     ctypedef enum LocalDescriptorOrientation:

freud/order.pxi

@@ -446,14 +446,17 @@ cdef class HexOrderParameter:
     :param k: symmetry of order parameter (:math:`k=6` is hexatic)
     :param n: number of neighbors (:math:`n=k` if :math:`n` not specified)
     :type rmax: float
-    :type k: unsigned int
+    :type k: float
     :type n: unsigned int
+
+    .. note:: While :math:`k` is a float, this is due to its use in calculations requiring floats. Passing in \
+    non-integer values will result in undefined behavior
     """
     cdef order.HexOrderParameter *thisptr
     cdef num_neigh
     cdef rmax
 
-    def __cinit__(self, rmax, k=int(6), n=int(0)):
+    def __cinit__(self, rmax, k=float(6.0), n=int(0)):
         self.thisptr = new order.HexOrderParameter(rmax, k, n)
         self.rmax = rmax
         self.num_neigh = (n if n else int(k))
@@ -553,7 +556,10 @@ cdef class HexOrderParameter:
         Get the symmetry of the order parameter
 
         :return: :math:`k`
-        :rtype: unsigned int
+        :rtype: float
+
+        .. note:: While :math:`k` is a float, this is due to its use in calculations requiring floats. Passing in \
+        non-integer values will result in undefined behavior
         """
         return self.getK()
 
@@ -562,9 +568,12 @@ cdef class HexOrderParameter:
         Get the symmetry of the order parameter
 
         :return: :math:`k`
-        :rtype: unsigned int
+        :rtype: float
+
+        .. note:: While :math:`k` is a float, this is due to its use in calculations requiring floats. Passing in \
+        non-integer values will result in undefined behavior
         """
-        cdef unsigned int k = self.thisptr.getK()
+        cdef float k = self.thisptr.getK()
         return k
 
 cdef class LocalDescriptors:

freud/order.py

@@ -6,19 +6,89 @@
 # Methods to compute order parameters
 #
 
-from ._freud import BondOrder
-from ._freud import CubaticOrderParameter
-from ._freud import HexOrderParameter
-from ._freud import TransOrderParameter
-from ._freud import LocalDescriptors
-from ._freud import Pairing2D
-from ._freud import AngularSeparation
-
-# everything below is spherical harmonic stuff
+# __all__ = ['HexOrderParameter']
+
+# not sure if broken
+from ._freud import BondOrder;
+from ._freud import CubaticOrderParameter;
+from ._freud import HexOrderParameter;
+from ._freud import TransOrderParameter;
+from ._freud import LocalDescriptors;
+from ._freud import Pairing2D;
+from ._freud import AngularSeparation;
+
+# everything below is sphericalharmonic stuff
 from ._freud import LocalQl
 from ._freud import LocalQlNear
 from ._freud import LocalWl
 from ._freud import LocalWlNear
 from ._freud import MatchEnv
 from ._freud import SolLiq
 from ._freud import SolLiqNear
+#import scipy.io as sio
+import numpy as np
+
+'''class LocalWl(LWl):
+
+    def __init__(self,box,rmax,l):
+        #Initialize those wigner3j values as mx1 numpy arra
+        # type of numpy should nd.type=float64??? check.
+        super(LocalWl,self).__init__(box,rmax,l)
+        #pywig3j = this.wigner3j(l);
+        if l < 22:
+            super(LocalWl,self).setWigner3j(self.getwigner(l))
+        else:
+            print('l too big, need sympy library')
+            super(LocalWl,self).setWigner3j(self.wigner3j(l))
+
+    #read of wigner3j coefficients from wigner3j.mat
+    def getwigner(self,l):
+        allwig = sio.loadmat('wigner3j.mat')
+        W_l = np.array(allwig['l'+str(l)][0],dtype=np.float64)
+        return W_l
+
+    #calculate wigner3j coefficients from sympy python library
+    def wigner3j(self,l):
+        from sympy.physics.wigner import wigner_3j
+        counter = 0
+        for u1 in range(2*l+1):
+            for u2 in range(max(l-u1,0),min(2*l+1,3*l-u1+1)):
+                counter += 1
+        W_l = np.zeros (counter,dtype=np.float64)
+
+        i = 0
+        for u1 in range(2*l+1):
+            for u2 in range(max(l-u1,0),min(2*l+1,3*l-u1+1)):
+                u3 = 3*l-u2-u1
+                W_l[i] = float(wigner_3j(l,l,l,u1-l,u2-l,u3-l))
+                i += 1
+        return W_l'''
+
+
+#How to set up wigner3j.mat file
+#The structure of this .mat file is a list of float numbers in the order of l
+# 1st row: l=2 wigner3j, 2nd row: l=4 wigner3j...
+#The coefficients are in the order of how the loop is written
+
+''' import scipy.io as sio
+ from sympy.physics.wigner import wigner_3j
+ import numpy as np
+
+ wigner={}
+ for i in range(10):
+     counter = 0
+     l = (i+1)*2
+     for u1 in range(2*l+1):
+         for u2 in range(max(l-u1,0),min(2*l+1,3*l-u1+1)):
+             counter += 1
+     W_l = np.zeros(counter,dtype=np.float64)
+
+     j = 0
+     for u1 in range(2*l+1):
+         for u2 in range(max(l-u1,0),min(2*l+1,3*l-u1+1)):
+             u3 = 3*l-u2-u1
+             W_l[j] = float(wigner_3j(l,l,l,u1-l,u2-l,u3-l))
+             j += 1
+     wigner['l'+str(l)]=W_l
+
+sio.savemat('wigner3j.mat',wigner)'''

tests/test_order_HexOrderParameter.py → tests/test_hexorder.py

@@ -12,17 +12,17 @@ def test_getK(self):
         box = freud.box.Box.square(boxlen)
 
         hop = freud.order.HexOrderParameter(rmax)
-        npt.assert_equal(hop.k, 6)
+        npt.assert_equal(hop.k, 6.0)
 
     def test_getK_pass(self):
         boxlen = 10
         N = 500
         rmax = 3
-        k = 3
+        k=3.0
 
         box = freud.box.Box.square(boxlen)
         hop = freud.order.HexOrderParameter(rmax, k)
-        npt.assert_equal(hop.k, 3)
+        npt.assert_equal(hop.k, 3.0)
 
     def test_getNP(self):
         boxlen = 10