From fd13756b868d1294f5fa68e9c7d1fb7f129af233 Mon Sep 17 00:00:00 2001 From: =?UTF-8?q?Marek=20Ne=C4=8Dada?= Date: Fri, 8 Mar 2019 11:10:01 +0000 Subject: [PATCH] ss vector packing and unpacking fixed and cython-wrapped. Former-commit-id: ab1002b4c7b6fa4b7c14520d689cd7cd030b84e0 --- qpms/qpms_c.pyx | 32 ++++++++++++++++++++++++++++++++ qpms/qpms_cdefs.pxd | 15 ++++++++++++++- qpms/scatsystem.c | 1 + tests/scatsys.py | 9 +++++++++ 4 files changed, 56 insertions(+), 1 deletion(-) diff --git a/qpms/qpms_c.pyx b/qpms/qpms_c.pyx index 1a1b7cd..dd498ce 100644 --- a/qpms/qpms_c.pyx +++ b/qpms/qpms_c.pyx @@ -1171,6 +1171,9 @@ cdef char *make_c_string(pythonstring): s[n] = 0 return s +def string_c2py(const char* cstring): + return cstring.decode('UTF-8') + cdef class FinitePointGroup: ''' Wrapper over the qpms_finite_group_t structure. @@ -1386,6 +1389,35 @@ cdef class ScatteringSystem: r.append(self.s[0].p[pi]) return r + property fecv_size: + def __get__(self): return self.s[0].fecv_size + property saecv_sizes: + def __get__(self): return [self.s[0].saecv_sizes[i] for i in range(self.s[0].sym[0].nirreps)] + property irrep_names: + def __get__(self): + return [string_c2py(self.s[0].sym[0].irreps[iri].name) if (self.s[0].sym[0].irreps[iri].name) else None + for iri in range(self.s[0].sym[0].nirreps)] + property nirreps: + def __get__(self): return self.s[0].sym[0].nirreps + + def pack_vector(self, vect, iri): + if len(vect) != self.fecv_size: raise ValueError("Length of a full vector has to be %d, not %d" % (self.fecv_size, len(vect))) + vect = np.array(vect, dtype=complex, copy=False, order='C') + cdef cdouble[::1] vect_view = vect; + cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty((self.saecv_sizes[iri],), dtype=complex) + cdef cdouble[::1] target_view = target_np + qpms_scatsys_irrep_pack_vector(&target_view[0], &vect_view[0], self.s, iri) + return target_np + def unpack_vector(self, packed, iri): + if len(packed) != self.saecv_sizes[iri]: raise ValueError("Length of %d. irrep-packed vector has to be %d, not %d" + % (iri, self.saecv_sizes, len(packed))) + packed = np.array(packed, dtype=complex, copy=False, order='C') + cdef cdouble[::1] packed_view = packed + cdef np.ndarray[np.complex_t, ndim=1] target_np = np.empty((self.fecv_size,), dtype=complex) + cdef cdouble[::1] target_view = target_np + qpms_scatsys_irrep_unpack_vector(&target_view[0], &packed_view[0], self.s, iri, 0) + return target_np + def tlm2uvswfi(t, l, m): ''' TODO doc And TODO this should rather be an ufunc. diff --git a/qpms/qpms_cdefs.pxd b/qpms/qpms_cdefs.pxd index e7f34fe..0d943f8 100644 --- a/qpms/qpms_cdefs.pxd +++ b/qpms/qpms_cdefs.pxd @@ -244,7 +244,10 @@ cdef extern from "scatsystem.h": qpms_ss_tmi_t tm_count qpms_particle_tid_t *p qpms_ss_pi_t p_count - # We shouldn't need more to construct a symmetric scatsystem + # We shouldn't need more to construct a symmetric scatsystem ^^^ + size_t fecv_size + size_t *saecv_sizes + const qpms_finite_group_t *sym qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) void qpms_scatsys_free(qpms_scatsys_t *s) qpms_errno_t qpms_scatsys_dump(qpms_scatsys_t *ss, char *path) #NI @@ -274,4 +277,14 @@ cdef extern from "scatsystem.h": qpms_errno_t qpms_load_scuff_tmatrix(const char *path, const qpms_vswf_set_spec_t *bspec, size_t *n, double **freqs, double **freqs_su, qpms_tmatrix_t **tmatrices_array, cdouble **tmdata) + cdouble *qpms_scatsys_irrep_pack_matrix(cdouble *target_packed, + const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) + cdouble *qpms_scatsys_irrep_unpack_matrix(cdouble *target_full, + const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add) + cdouble *qpms_scatsys_irrep_pack_vector(cdouble *target_packed, + const cdouble *orig_full, const qpms_scatsys_t *ss, qpms_iri_t iri) + cdouble *qpms_scatsys_irrep_unpack_vector(cdouble *target_full, + const cdouble *orig_packed, const qpms_scatsys_t *ss, qpms_iri_t iri, bint add) + + diff --git a/qpms/scatsystem.c b/qpms/scatsystem.c index b075605..fecc385 100644 --- a/qpms/scatsystem.c +++ b/qpms/scatsystem.c @@ -405,6 +405,7 @@ static void add_orbit_type(qpms_scatsys_t *ss, const qpms_ss_orbit_type_t *ot_cu "expected %d = %d * %d, got %d.", ot_new->size * bspecn, ot_new->size, bspecn, bs_cumsum); ot_new->instance_count = 0; + ss->orbit_type_count++; } diff --git a/tests/scatsys.py b/tests/scatsys.py index 3367200..f42d6e9 100644 --- a/tests/scatsys.py +++ b/tests/scatsys.py @@ -11,3 +11,12 @@ p1 = Particle((1,2,),t1) p2 = Particle((1,2,3),t1) p3 = Particle((0.1,2),t2) ss = ScatteringSystem([p1, p2, p3], sym) + +#print(ss.fecv_size, ss.saecv_sizes, ss.nirreps, ss.nirrep_names) + +fullvector = np.random.rand(ss.fecv_size) + 1j*np.random.rand(ss.fecv_size) +packedvectors = [(iri, ss.pack_vector(fullvector, iri)) for iri in range(ss.nirreps)] +unpackedvectors = np.array([ss.unpack_vector(v[1], v[0]) for v in packedvectors]) +rec_fullvector = np.sum(unpackedvectors, axis=0) +thediff = np.amax(abs(rec_fullvector-fullvector)) +assert(thediff < 1e-14)