ss vector packing and unpacking fixed and cython-wrapped.

Former-commit-id: ab1002b4c7b6fa4b7c14520d689cd7cd030b84e0

qpms/qpms_c.pyx

@@ -1171,6 +1171,9 @@ cdef char *make_c_string(pythonstring):
     s[n] = <char>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.

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)
+
+
 

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++;
 }
 
 

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)