PointGroup cython class, fix many bugs (quaternion ambiguity).

Former-commit-id: 16835b9f03c4fbb28c414ce9844cabb30aabfb00

qpms/pointgroups.c

@@ -7,12 +7,13 @@
   if ((a) > (b)) return 1;\
 }
 
-int qpms_pg_irot3_cmp(const qpms_irot3_t *a, const qpms_irot3_t *b) {
+int qpms_pg_irot3_cmp(const qpms_irot3_t *p1, const qpms_irot3_t *p2) {
-  PAIRCMP(a->det, b->det);
+  PAIRCMP(p1->det, p2->det);
-  PAIRCMP(creal(a->rot.a), creal(b->rot.a));
+  const qpms_quat_t r1 = qpms_quat_standardise(p1->rot), r2 = qpms_quat_standardise(p2->rot);
-  PAIRCMP(cimag(a->rot.a), cimag(b->rot.a));
+  PAIRCMP(creal(r1.a), creal(r2.a));
-  PAIRCMP(creal(a->rot.b), creal(b->rot.b));
+  PAIRCMP(cimag(r1.a), cimag(r2.a));
-  PAIRCMP(cimag(a->rot.b), cimag(b->rot.b));
+  PAIRCMP(creal(r1.b), creal(r2.b));
+  PAIRCMP(cimag(r1.b), cimag(r2.b));
   return 0;
 }
 
@@ -38,16 +39,15 @@ int qpms_pg_irot3_approx_cmp_v(const void *av, const void *bv) {
 qpms_irot3_t *qpms_pg_canonical_elems(qpms_irot3_t *target,
     qpms_pointgroup_class cls, const qpms_gmi_t then) {
   QPMS_UNTESTED;
-  qpms_gmi_t order = qpms_pg_order(cls, then);
+  const qpms_gmi_t order = qpms_pg_order(cls, then);
   QPMS_ENSURE(order, "Cannot generate an infinite group!");
-  if (!target) QPMS_CRASHING_MALLOC(target, order * sizeof(qpms_irot3_t));
+  if (!target) QPMS_CRASHING_MALLOC(target, (1+order) * sizeof(qpms_irot3_t));
   target[0] = QPMS_IROT3_IDENTITY;
   qpms_gmi_t ngen = qpms_pg_genset_size(cls, then);
   qpms_irot3_t gens[ngen]; 
   (void) qpms_pg_genset(cls, then, gens);
 
   // Let's try it with a binary search tree, as an exercise :)
-  qpms_irot3_t tree[order];
   void *root = NULL;
   // put the starting element (identity) to the tree
   (void) tsearch((void *) target, &root, qpms_pg_irot3_approx_cmp_v);
@@ -62,7 +62,7 @@ qpms_irot3_t *qpms_pg_canonical_elems(qpms_irot3_t *target,
 
   while (si >= 0) { // DFS
     if (gistack[si] < ngen) { // are there generators left at this level? If so, try another elem
-      if (n >= order) QPMS_WTF; // TODO some error message
+      if (n > order) QPMS_WTF; // TODO some error message
       target[n] = qpms_irot3_mult(gens[gistack[si]], target[srcstack[si]]);
       if (tfind((void *) &(target[n]), &root, qpms_pg_irot3_approx_cmp_v))
         // elem found, try it with another gen in the next iteration
@@ -85,7 +85,7 @@ qpms_irot3_t *qpms_pg_canonical_elems(qpms_irot3_t *target,
       "(assumed group order = %d, got %d; qpms_pg_quat_cmp_atol = %g)",
       order, n, qpms_pg_quat_cmp_atol);
 
-  while(root) tdelete(root, &root, qpms_pg_irot3_approx_cmp_v); // I hope this is the correct way.
+  while(root) tdelete(*(qpms_irot3_t **)root, &root, qpms_pg_irot3_approx_cmp_v); // I hope this is the correct way.
 
   return target;
 }

qpms/qpms_c.pyx

@@ -32,6 +32,29 @@ class BesselType(enum.IntEnum):
     HANKEL_PLUS = QPMS_HANKEL_PLUS
     HANKEL_MINUS = QPMS_HANKEL_MINUS
 
+class PointGroupClass(enum.IntEnum):
+    CN    = QPMS_PGS_CN
+    S2N   = QPMS_PGS_S2N
+    CNH   = QPMS_PGS_CNH
+    CNV   = QPMS_PGS_CNV
+    DN    = QPMS_PGS_DN
+    DND   = QPMS_PGS_DND
+    DNH   = QPMS_PGS_DNH
+    T     = QPMS_PGS_T
+    TD    = QPMS_PGS_TD
+    TH    = QPMS_PGS_TH
+    O     = QPMS_PGS_O
+    OH    = QPMS_PGS_OH
+    I     = QPMS_PGS_I
+    IH    = QPMS_PGS_IH
+    CINF  = QPMS_PGS_CINF
+    CINFH = QPMS_PGS_CINFH
+    CINFV = QPMS_PGS_CINFV
+    DINF  = QPMS_PGS_DINF
+    DINFH = QPMS_PGS_DINFH
+    SO3   = QPMS_PGS_SO3
+    O3    = QPMS_PGS_O3
+
 class VSWFNorm(enum.IntEnum):
     # TODO try to make this an enum.IntFlag if supported
     # TODO add the other flags from qpms_normalisation_t as well
@@ -1022,6 +1045,9 @@ cdef class IRot3:
         else:
             raise ValueError('Unsupported constructor arguments')
 
+    cdef void cset(self, qpms_irot3_t qd):
+        self.qd = qd
+
     def copy(self):
         res = IRot3(CQuat(1,0,0,0),1)
         res.qd = self.qd
@@ -1329,6 +1355,49 @@ cdef char *make_c_string(pythonstring):
 def string_c2py(const char* cstring):
     return cstring.decode('UTF-8')
 
+cdef class PointGroup:
+    cdef readonly qpms_pointgroup_t G
+
+    def __init__(self, cls, qpms_gmi_t n = 0, IRot3 orientation = IRot3()):
+        cls = PointGroupClass(cls)
+        self.G.c = cls
+        if n <= 0 and qpms_pg_is_finite_axial(cls):
+            raise ValueError("For finite axial groups, n argument must be positive")
+        self.G.n = n
+        self.G.orientation = orientation.qd
+
+    def __len__(self):
+        return qpms_pg_order(self.G.c, self.G.n);
+
+    def __le__(PointGroup self, PointGroup other):
+        if qpms_pg_is_subgroup(self.G, other.G):
+            return True
+        else:
+            return False
+    def __ge__(PointGroup self, PointGroup other):
+        if qpms_pg_is_subgroup(other.G, self.G):
+            return True
+        else:
+            return False
+    def __lt__(PointGroup self, PointGroup other):
+        return qpms_pg_is_subgroup(self.G, other.G) and not qpms_pg_is_subgroup(other.G, self.G)
+    def __eq__(PointGroup self, PointGroup other):
+        return qpms_pg_is_subgroup(self.G, other.G) and qpms_pg_is_subgroup(other.G, self.G)
+    def __gt__(PointGroup self, PointGroup other):
+        return not qpms_pg_is_subgroup(self.G, other.G) and qpms_pg_is_subgroup(other.G, self.G)
+
+    def elems(self):
+        els = list()
+        cdef qpms_irot3_t *arr
+        arr = qpms_pg_elems(NULL, self.G)
+        cdef IRot3 q
+        for i in range(len(self)):
+            q = IRot3()
+            q.cset(arr[i])
+            els.append(q)
+        free(arr)
+        return els
+
 cdef class FinitePointGroup:
     '''
     Wrapper over the qpms_finite_group_t structure.

qpms/qpms_cdefs.pxd

@@ -103,6 +103,32 @@ cdef extern from "qpms_types.h":
         qpms_vswf_set_spec_t *spec
         cdouble *m
         bint owns_m # FIXME in fact bool
+    ctypedef enum qpms_pointgroup_class:
+        QPMS_PGS_CN
+        QPMS_PGS_S2N
+        QPMS_PGS_CNH
+        QPMS_PGS_CNV
+        QPMS_PGS_DN
+        QPMS_PGS_DND
+        QPMS_PGS_DNH
+        QPMS_PGS_T
+        QPMS_PGS_TD
+        QPMS_PGS_TH
+        QPMS_PGS_O
+        QPMS_PGS_OH
+        QPMS_PGS_I
+        QPMS_PGS_IH
+        QPMS_PGS_CINF
+        QPMS_PGS_CINFH
+        QPMS_PGS_CINFV
+        QPMS_PGS_DINF
+        QPMS_PGS_DINFH
+        QPMS_PGS_SO3
+        QPMS_PGS_O3
+    struct qpms_pointgroup_t:
+        qpms_pointgroup_class c
+        qpms_gmi_t n
+        qpms_irot3_t orientation
     # maybe more if needed
 
 cdef extern from "qpms_error.h":
@@ -349,6 +375,20 @@ cdef extern from "tmatrices.h":
     double qpms_permittivity_interpolator_omega_min(const qpms_permittivity_interpolator_t *interp)
     void qpms_permittivity_interpolator_free(qpms_permittivity_interpolator_t *interp)
 
+cdef extern from "pointgroups.h":
+    bint qpms_pg_is_finite_axial(qpms_pointgroup_class cls)
+    double qpms_pg_quat_cmp_atol
+    int qpms_pg_irot3_cmp(const qpms_irot3_t *, const qpms_irot3_t *);
+    int qpms_pg_irot3_cmp_v(const void *, const void *);
+    int qpms_pg_irot3_approx_cmp(const qpms_irot3_t *a, const qpms_irot3_t *b, double atol)
+    int qpms_pg_irot3_approx_cmp_v(const void *a, const void *b)
+
+    qpms_gmi_t qpms_pg_order(qpms_pointgroup_class cls, qpms_gmi_t n)
+    qpms_irot3_t *qpms_pg_canonical_elems( qpms_irot3_t *target, qpms_pointgroup_class cls, qpms_gmi_t n)
+    qpms_gmi_t qpms_pg_genset_size(qpms_pointgroup_class cls, qpms_gmi_t n)
+    qpms_gmi_t qpms_pg_genset(qpms_pointgroup_class cls, qpms_gmi_t n, qpms_irot3_t *gen)
+    qpms_irot3_t *qpms_pg_elems(qpms_irot3_t *target, qpms_pointgroup_t g)
+    bint qpms_pg_is_subgroup(qpms_pointgroup_t a, qpms_pointgroup_t b);
 
 cdef extern from "scatsystem.h":
     void qpms_scatsystem_set_nthreads(long n)
@@ -415,3 +455,4 @@ cdef extern from "scatsystem.h":
             int *target_piv, const qpms_scatsys_t *ss, qpms_iri_t iri)
     cdouble *qpms_scatsys_scatter_solve(cdouble *target_f, const cdouble *a_inc, qpms_ss_LU ludata)
 
+

qpms/quaternions.h

@@ -97,6 +97,33 @@ static inline bool qpms_quat_isclose(const qpms_quat_t p, const qpms_quat_t q, d
 	return qpms_quat_norm(qpms_quat_sub(p,q)) <= atol;
 }
 
+/// "Standardises" a quaternion to have the largest component "positive".
+/**
+ * This is to remove the ambiguity stemming from the double cover of SO(3).
+ */
+static inline qpms_quat_t qpms_quat_standardise(qpms_quat_t p) {
+	double maxabs = 0;
+	int maxi = 0;
+	const double *arr = (double *) &(p.a);
+	for(int i = 0; i < 4; ++i)
+		if (fabs(arr[i]) > maxabs) {
+			maxi = i;
+			maxabs = fabs(arr[i]);
+		}
+	if(arr[maxi] < 0) {
+		p.a = -p.a;
+		p.b = -p.b;
+	}
+	return p;
+}		
+
+/// Test approximate equality of "standardised" quaternions, so that \f$-q\f$ is considered equal to \f$q\f$.
+static inline bool qpms_quat_isclose2(const qpms_quat_t p, const qpms_quat_t q, double atol) {
+	return qpms_quat_norm(qpms_quat_sub(
+				qpms_quat_standardise(p),
+				qpms_quat_standardise(q))) <= atol;
+}
+
 /// Norm of the quaternion imaginary (vector) part.
 static inline double qpms_quat_imnorm(const qpms_quat_t q) {
 	const double z = cimag(q.a), x = cimag(q.b), y = creal(q.b);
@@ -224,7 +251,7 @@ static inline qpms_irot3_t qpms_irot3_pow(const qpms_irot3_t p, int n) {
 
 /// Test approximate equality of irot3.
 static inline bool qpms_irot3_isclose(const qpms_irot3_t p, const qpms_irot3_t q, double atol) {
-	return qpms_quat_isclose(p.rot, q.rot, atol) && p.det == q.det;
+	return qpms_quat_isclose2(p.rot, q.rot, atol) && p.det == q.det;
 }
 
 /// Apply an improper rotation onto a 3d cartesian vector.
@@ -261,7 +288,7 @@ static inline qpms_quat_t qpms_quat_zrot_angle(double angle) {
 
 /// versor representing rotation \f$ C_N^k \f$, i.e. of angle \f$ 2\pi k / N\f$ around z axis.
 static inline qpms_quat_t qpms_quat_zrot_Nk(double N, double k) {
-	return qpms_quat_zrot_angle(M_PI * k / N);
+	return qpms_quat_zrot_angle(2 * M_PI * k / N);
 }
 
 /// Rotation around z-axis.
@@ -272,7 +299,7 @@ static inline qpms_irot3_t qpms_irot3_zrot_angle(double angle) {
 
 /// Rotation \f$ C_N^k \f$, i.e. of angle \f$ 2\pi k / N\f$ around z axis.
 static inline qpms_irot3_t qpms_irot3_zrot_Nk(double N, double k) {
-	return qpms_irot3_zrot_angle(M_PI * k / N);
+	return qpms_irot3_zrot_angle(2 * M_PI * k / N);
 }
 
 #endif //QPMS_WIGNER_H

setup.py

@@ -75,8 +75,9 @@ qpms_c = Extension('qpms_c',
             'qpms/error.c',
             'qpms/bessel.c',
             'qpms/own_zgemm.c',
+            'qpms/pointgroups.c',
         ],
-        extra_compile_args=['-std=c99','-ggdb', '-O3',
+        extra_compile_args=['-std=c99','-ggdb', '-O0',
             '-DQPMS_COMPILE_PYTHON_EXTENSIONS', # this is required
             #'-DQPMS_USE_OMP',
             '-DQPMS_SCATSYSTEM_USE_OWN_BLAS',

tests/pointgrouptest.py

@@ -0,0 +1,27 @@
+#!/usr/bin/env python3
+from qpms import PointGroup, CQuat, IRot3, PointGroupClass as PGC
+
+try:
+    PointGroup(PGC.CN)
+except ValueError:
+    print("OK")
+
+try:
+    PointGroup(66)
+except ValueError:
+    print("OK")
+
+C1 = PointGroup(PGC.CN, 1)
+print(len(C1))
+print(C1.elems())
+
+
+C2 = PointGroup(PGC.CN, 2)
+print(len(C2))
+print(C2.elems())
+
+D2H = PointGroup(PGC.DNH, 2)
+print(len(D2H))
+print(D2H.elems())
+
+exit(0)