From 1b63d7515363a1f536be54f4f028ff6ee5d71086 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Wed, 24 Jul 2019 18:26:12 +0300 Subject: [PATCH] PointGroup cython class, fix many bugs (quaternion ambiguity). Former-commit-id: 16835b9f03c4fbb28c414ce9844cabb30aabfb00 --- qpms/pointgroups.c | 22 ++++++------- qpms/qpms_c.pyx | 69 +++++++++++++++++++++++++++++++++++++++++ qpms/qpms_cdefs.pxd | 41 ++++++++++++++++++++++++ qpms/quaternions.h | 33 ++++++++++++++++++-- setup.py | 3 +- tests/pointgrouptest.py | 27 ++++++++++++++++ 6 files changed, 180 insertions(+), 15 deletions(-) create mode 100644 tests/pointgrouptest.py diff --git a/qpms/pointgroups.c b/qpms/pointgroups.c index c76fe95..063768b 100644 --- a/qpms/pointgroups.c +++ b/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) { - PAIRCMP(a->det, b->det); - PAIRCMP(creal(a->rot.a), creal(b->rot.a)); - PAIRCMP(cimag(a->rot.a), cimag(b->rot.a)); - PAIRCMP(creal(a->rot.b), creal(b->rot.b)); - PAIRCMP(cimag(a->rot.b), cimag(b->rot.b)); +int qpms_pg_irot3_cmp(const qpms_irot3_t *p1, const qpms_irot3_t *p2) { + PAIRCMP(p1->det, p2->det); + const qpms_quat_t r1 = qpms_quat_standardise(p1->rot), r2 = qpms_quat_standardise(p2->rot); + PAIRCMP(creal(r1.a), creal(r2.a)); + PAIRCMP(cimag(r1.a), cimag(r2.a)); + 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; } diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index a56a85a..a8e5960 100644 --- a/qpms/qpms_c.pyx +++ b/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. diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index e9d5b70..682f92f 100644 --- a/qpms/qpms_cdefs.pxd +++ b/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) + diff --git a/qpms/quaternions.h b/qpms/quaternions.h index 22c8643..fb1e282 100644 --- a/qpms/quaternions.h +++ b/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 diff --git a/setup.py b/setup.py index 8578bec..40b0a87 100755 --- a/setup.py +++ b/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', diff --git a/tests/pointgrouptest.py b/tests/pointgrouptest.py new file mode 100644 index 0000000..a14eadf --- /dev/null +++ b/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)