Towards a new scattering system (scatsystem.c compiles, but with some

warnings).


Former-commit-id: 8e873fa9eeed41d5384aac4844fdebf768063fe8

qpms/groups.h

@@ -19,7 +19,7 @@ typedef const char * qpms_permutation_t;
 /// To be used only in qpms_finite_group_t
 struct qpms_finite_group_irrep_t {
 	int dim; ///< Irrep dimension.
-	char name[]; ///< Irrep label.
+	char *name; ///< Irrep label.
 	/// Irrep matrix data.
 	/** The r-th row, c-th column of the representation of the i'th element is retrieved as
 	 *  m[i * dim * dim + r * dim + c]
@@ -49,12 +49,12 @@ typedef struct qpms_finite_group_t {
 	qpms_gmi_t *invi; ///< Group elem inverse indices.
 	qpms_gmi_t *gens; ///< A canonical set of group generators.
 	int ngens; ///< Number of the generators in gens;
-	qpms_permutation_t permrep[]; ///< Permutation representations of the elements.
+	qpms_permutation_t *permrep; ///< Permutation representations of the elements.
 	char **elemlabels; ///< Optional human readable labels for the group elements.
 	int permrep_nelem; ///< Number of the elements over which the permutation representation acts.
-	struct qpms_irot3_t rep3d[]; ///< The quaternion representation of a 3D point group (if applicable).
+	struct qpms_irot3_t *rep3d; ///< The quaternion representation of a 3D point group (if applicable).
 	int nirreps; ///< How many irreps does the group have
-	struct qpms_finite_group_irrep_t irreps[]; ///< Irreducible representations of the group.
+	struct qpms_finite_group_irrep_t *irreps; ///< Irreducible representations of the group.
 } qpms_finite_group_t;
 
 

qpms/qpms_cdefs.pxd

@@ -185,5 +185,18 @@ cdef extern from "translations.h":
             char *phi_data, np.npy_intp *phi_shape, np.npy_intp *phi_strides,
             char *r_ge_d_data, np.npy_intp *phi_shape, np.npy_intp *phi_strides) nogil
 
+cdef extern from "scatsystem.h":
+    struct qpms_tmatrix_t:
+        qpms_vswf_set_spec_t *spec
+        cdouble *m
+        int owns_m # FIXME in fact bool
+    struct qpms_particle_t:
+        cart3_t pos
+        const qpms_tmatrix_t *tmatrix
+    struct qpms_tmatrix_interpolator_t:
+        const qpms_vswf_set_spec_t *bspec
+    struct qpms_scatsys_t:
+        pass # TODO
+
 
 

qpms/scatsystem.c

@@ -2,8 +2,19 @@
 #include <cblas.h>
 #include "scatsystem.h"
 #include "indexing.h"
+#include "vswf.h"
+#include "groups.h"
+#include "symmetries.h"
 #include <gsl/gsl_spline.h>
+#include <assert.h>
+#include <unistd.h>
+#include "vectors.h"
+#include "wigner.h"
 
+
+#define QPMS_SCATSYS_LEN_RTOL 1e-13
+#define QPMS_SCATSYS_TMATRIX_ATOL 1e-14
+#define QPMS_SCATSYS_TMATRIX_RTOL 1e-12
 qpms_tmatrix_t *qpms_tmatrix_init(const qpms_vswf_set_spec_t *bspec) {
   qpms_tmatrix_t *t = malloc(sizeof(qpms_tmatrix_t));
   if (!t) abort();
@@ -30,6 +41,28 @@ void qpms_tmatrix_free(qpms_tmatrix_t *t){
   free(t);
 }
 
+qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
+                qpms_tmatrix_t *T, 
+                const complex double *M 
+                )
+{
+  //qpms_tmatrix_t *t = qpms_tmatrix_init(T->spec);
+  const size_t n = T->spec->n;
+  complex double tmp[n][n];
+  // tmp = M T
+  const complex double one = 1, zero = 0;
+  cblas_zgemm(CblasRowMajor,
+      CblasNoTrans,
+      CblasNoTrans,
+      n, n, n, &one, M, n, T->m, n, &zero, tmp, n);
+  // t->m = tmp M* = M T M*
+  cblas_zgemm(CblasRowMajor,
+      CblasNoTrans,
+      CblasConjTrans,
+      n, n, n, &one, tmp, n, M, n, &zero, T->m, n);
+  return T;
+}
+
 qpms_tmatrix_t *qpms_tmatrix_apply_symop(
                 const qpms_tmatrix_t *T, 
                 const complex double *M 
@@ -129,7 +162,7 @@ qpms_tmatrix_t *qpms_tmatrix_symmetrise_C_N_inplace(qpms_tmatrix_t *T, int N) {
   return T;
 }
 
-bool qpms_tmatrix_isnear(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B,
+bool qpms_tmatrix_isclose(const qpms_tmatrix_t *A, const qpms_tmatrix_t *B,
     const double rtol, const double atol)
 {
   if (!qpms_vswf_set_spec_isidentical(A->spec, B->spec))
@@ -218,3 +251,146 @@ qpms_tmatrix_t *qpms_tmatrix_interpolator_eval(const qpms_tmatrix_interpolator_t
   return t;
 }
 
+
+
+// ------------ Stupid implementation of qpms_scatsys_apply_symmetry() -------------
+
+#define MIN(x,y) (((x)<(y))?(x):(y))
+#define MAX(x,y) (((x)>(y))?(x):(y))
+
+qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const qpms_finite_group_t *sym) {
+  // TODO check data sanity
+
+  // First, determine the rough radius of the array
+  double lenscale = 0;
+  {
+    double minx = +INFINITY, miny = +INFINITY, minz = +INFINITY;
+    double maxx = -INFINITY, maxy = -INFINITY, maxz = -INFINITY;
+    for (ssize_t i = 0; i < orig->p_count; ++i) {
+      minx = MIN(minx,orig->p[i].pos.x);
+      miny = MIN(miny,orig->p[i].pos.y);
+      minz = MIN(minz,orig->p[i].pos.z);
+      maxx = MAX(maxx,orig->p[i].pos.x);
+      maxy = MAX(maxy,orig->p[i].pos.y);
+      maxz = MAX(maxz,orig->p[i].pos.z);
+    }
+    lenscale = ((maxx-minx)+(maxy-miny)+(maxz-minz)) / 3;
+  }
+
+  // Second, check that there are no duplicit positions in the input system.
+  for (ssize_t i = 0; i < orig->p_count; ++i)
+    for (ssize_t j = 0; j < i; ++j)
+      assert(!cart3_isclose(orig->p[i].pos, orig->p[j].pos, 0, QPMS_SCATSYS_LEN_RTOL * lenscale));
+
+  // Allocate T-matrix and particle arrays
+  qpms_scatsys_t *ss = malloc(sizeof(qpms_scatsys_t));
+  ss->lenscale = lenscale;
+
+  ss->tm_capacity = sym->order * orig->tm_count;
+  ss->tm = malloc(ss->tm_capacity * sizeof(qpms_tmatrix_t *));
+
+  ss->p_capacity = sym->order * orig->p_count;
+  ss->p = malloc(ss->p_capacity * sizeof(qpms_particle_tid_t));
+
+  // Copy T-matrices; checking for duplicitie
+
+  ptrdiff_t tm_dupl_remap[ss->tm_capacity]; // Auxilliary array to label remapping the indices after ignoring t-matrix duplicities
+  ss->tm_count = 0;
+  for (ptrdiff_t i = 0; i < orig->tm_count; ++i) {
+    bool found_dupl = false;
+    ptrdiff_t j;
+    for (j = 0; j < ss->tm_count; ++j) 
+      if (qpms_tmatrix_isclose(orig->tm[i], ss->tm[j], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL)) {
+        break;
+      }
+    if (j == ss->tm_count) { // duplicity not found, copy the t-matrix
+      ss->tm[i] = qpms_tmatrix_copy(orig->tm[i]);
+      ++(ss->tm_count);
+    } 
+    tm_dupl_remap[i] = j;
+  }
+
+  // Copy particles, remapping the t-matrix indices
+  for (ptrdiff_t i = 0; i < orig->p_count; ++(i)) {
+    ss->p[i] = orig->p[i];
+    ss->p[i].tmatrix_id = tm_dupl_remap[ss->p[i].tmatrix_id];
+  }
+  ss->p_count = orig->p_count;
+
+  // allocate t-matrix symmetry map
+  ss->tm_sym_map = malloc(sizeof(ptrdiff_t) * sym->order * sym->order * ss->tm_count);
+
+  // Extend the T-matrices list by the symmetry operations
+  for (ptrdiff_t tmi = 0; tmi < ss->tm_count; ++tmi) 
+    for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi){
+      const size_t d = ss->tm[tmi]->spec->n;
+      complex double M[d][d]; // transformation matrix
+      qpms_irot3_uvswfi_dense(M[0], ss->tm[tmi]->spec, sym->rep3d[gmi]);
+      qpms_tmatrix_t *transformed = qpms_tmatrix_apply_symop(ss->tm[tmi], M[0]);
+      ptrdiff_t tmj;
+      for (tmj = 0; tmj < ss->tm_count; ++tmj)
+        if (qpms_tmatrix_isclose(transformed, ss->tm[tmj], QPMS_SCATSYS_TMATRIX_RTOL, QPMS_SCATSYS_TMATRIX_ATOL))
+          break;
+      if (tmj < ss->tm_count) { // HIT, transformed T-matrix already exists
+        qpms_tmatrix_free(transformed);
+      } else { // MISS, save the matrix and increment the count
+        ss->tm[ss->tm_count] = transformed;
+        ++(ss->tm_count);
+      }
+      ss->tm_sym_map[gmi + tmi * sym->order] = tmj; // In any case, tmj now indexes the correct transformed matrix
+    }
+  // Possibly free some space using the new ss->tm_count instead of (old) ss->tm_count*sym->order
+  ss->tm_sym_map = realloc(ss->tm_sym_map, sizeof(ptrdiff_t) * sym->order * ss->tm_count);
+  // tm could be realloc'd as well, but those are just pointers, not likely many.
+  
+  // allocate particle symmetry map
+  ss->p_sym_map = malloc(sizeof(ptrdiff_t) * sym->order * sym->order * ss->p_count);
+  // Extend the particle list by the symmetry operations, check that particles mapped by symmetry ops on themselves
+  // have the correct symmetry
+  // TODO this could be sped up to O(npart * log(npart)); let's see later if it is needed.
+  for (ptrdiff_t pi = 0; pi < ss->p_count; ++pi)
+    for (qpms_gmi_t gmi = 0; gmi < sym->order; ++gmi) {
+      cart3_t newpoint = qpms_irot3_apply_cart3(sym->rep3d[gmi], ss->p[pi].pos);
+      ptrdiff_t new_tmi = ss->tm_sym_map[gmi + ss->p[pi].tmatrix_id * sym->order]; // transformed t-matrix index
+      ptrdiff_t pj;
+      for (pj = 0; pj < ss->p_count; ++pj)
+        if (cart3_isclose(newpoint, ss->p[pj].pos, 0, ss->lenscale * QPMS_SCATSYS_LEN_RTOL)) {
+          if (new_tmi != ss->p[pj].tmatrix_id)
+            abort(); // The particle is mapped to another (or itself) without having the required symmetry.  TODO give some error message.
+          break;
+        }
+      if (pj < ss->p_count) { // HIT, the particle is transformed to an "existing" one.
+        ;
+      } else { // MISS, the symmetry transforms the particle to a new location, so add it.
+        qpms_particle_tid_t newparticle = {newpoint, new_tmi};
+        ss->p[ss->p_count] = newparticle;
+        ++(ss->p_count);
+      }
+      ss->p_sym_map[gmi + pi * sym->order] = pi;
+    }
+  // Possibly free some space using the new ss->p_count instead of (old) ss->p_count*sym->order
+  ss->p_sym_map = realloc(ss->p_sym_map, sizeof(ptrdiff_t) * sym->order * ss->p_count);
+  ss->p = realloc(ss->p, sizeof(qpms_particle_tid_t) * ss->p_count); ss->p_capacity = ss->p_count;
+
+  ss->sym = sym;
+  return ss;
+}
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+
+  
+
+

qpms/scatsystem.h

@@ -52,7 +52,7 @@ typedef struct qpms_particle_t {
  *  This function actually checks for identical vswf specs.
  * TODO define constants with "default" atol, rtol for this function.
  */
-bool qpms_tmatrix_isnear(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
+bool qpms_tmatrix_isclose(const qpms_tmatrix_t *T1, const qpms_tmatrix_t *T2,
 		const double rtol, const double atol);
 
 /// Creates a T-matrix from another matrix and a symmetry operation.
@@ -66,6 +66,17 @@ qpms_tmatrix_t *qpms_tmatrix_apply_symop(
 		const complex double *M ///< the symmetry op matrix in row-major format
 		);
 
+/// Applies a symmetry operation onto a T-matrix, rewriting the original T-matrix data.
+/** The symmetry operation is expected to be a unitary (square) 
+ *  matrix \a M with the same
+ *  VSWF basis spec as the T-matrix (i.e. \a t->spec). The new T-matrix will then
+ *  correspond to CHECKME \f[ T' = MTM^\dagger \f] 
+ */
+qpms_tmatrix_t *qpms_tmatrix_apply_symop_inplace(
+		qpms_tmatrix_t *T, ///< the original T-matrix
+		const complex double *M ///< the symmetry op matrix in row-major format
+		);
+
 /// Symmetrizes a T-matrix with an involution symmetry operation.
 /** The symmetry operation is expected to be a unitary (square) 
  *  matrix \a M with the same
@@ -154,7 +165,7 @@ extern const gsl_interp_type * gsl_interp_akima_periodic;
 extern const gsl_interp_type * gsl_interp_steffen;
 */
 
-// struct gsl_interp_accel; // use if lookup shows to be too slow
+// struct gsl_interp_accel; // use if lookup proves to be too slow
 typedef struct qpms_tmatrix_interpolator_t {
 	const qpms_vswf_set_spec_t *bspec;
 	//bool owns_bspec;
@@ -177,8 +188,42 @@ qpms_tmatrix_interpolator_t *qpms_tmatrix_interpolator_create(size_t n, ///< Num
 	       //, bool copy_bspec ///< if true, copies its own copy of basis spec from the first T-matrix.
        	       /*, ...? */);
 
+
+/// A particle, defined by its T-matrix INDEX and position, to be used in qpms_scatsys_t.
+typedef struct qpms_particle_tid_t {
+	// Does it make sense to ever use other than cartesian coords for this?
+	cart3_t pos; ///< Particle position in cartesian coordinates.
+	ptrdiff_t tmatrix_id; ///< T-matrix index
+} qpms_particle_tid_t;
+
+struct qpms_finite_group_t;
+
 typedef struct qpms_scatsys_t {
-	;
+	// TODO does bspec belong here?
+	qpms_tmatrix_t **tm; ///< T-matrices in the system
+	ptrdiff_t tm_count; ///< Number of all different T-matrices
+	ptrdiff_t tm_capacity; ///< Capacity of tm[].
+	qpms_particle_tid_t *p; ///< Particles.
+	ptrdiff_t p_count; ///< Size of particles array.
+	ptrdiff_t p_capacity; ///< Capacity of p[].
+
+	//TODO the index types do not need to be so big.
+	struct qpms_finite_group_t *sym; ///< Symmetry group of the array
+	ptrdiff_t *p_sym_map; ///< Which particles map onto which by the symmetry ops.
+	///< p_sym_map[idi + pi * sym->order] gives the index of pi-th particle under the idi'th sym op.
+	ptrdiff_t *tm_sym_map; ///< Which t-matrices map onto which by the symmetry ops. Lookup by tm_sum_map[idi + tmi * sym->order].
+
+	// TODO shifted origin of the symmetry group etc.
+	// TODO some indices for fast operations here.
+	// private
+	double lenscale; // radius of the array, used as a relative tolerance measure
 } qpms_scatsys_t;
 
+/// Creates a new scatsys by applying a symmetry group, copying particles if needed.
+/** In fact, it copies everything except the vswf set specs, so keep them alive until scatsys is destroyed.
+ */
+qpms_scatsys_t *qpms_scatsys_apply_symmetry(const qpms_scatsys_t *orig, const struct qpms_finite_group_t *sym);
+/// Destroys the result of qpms_scatsys_apply_symmetry.
+void qpms_scatsys_free(qpms_scatsys_t *s);
+
 #endif //QPMS_SCATSYSTEM_H

qpms/vectors.h

@@ -122,6 +122,14 @@ static inline cart3_t cart3_scale(const double c, const cart3_t v) {
 	return res;
 }
 
+static inline double cart3_dist(const cart3_t a, const cart3_t b) {
+	return cart3norm(cart3_substract(a,b));
+}
+
+static inline bool cart3_isclose(const cart3_t a, const cart3_t b, double rtol, double atol) {
+	return cart3_dist(a,b) <= atol + rtol * (cart3norm(b) + cart3norm(a)) * .5;
+}
+
 static inline ccart3_t ccart3_scale(const complex  double c, const ccart3_t v) {
 	ccart3_t res = {c * v.x, c * v.y, c * v.z};
 	return res;
@@ -161,6 +169,7 @@ static inline double cart3_dot(const cart3_t a, const cart3_t b) {
 	return a.x * b.x + a.y * b.y + a.z * b.z;
 }
 
+
 // equivalent to sph_loccart2cart in qpms_p.py
 static inline ccart3_t csphvec2ccart(const csphvec_t sphvec, const sph_t at) {
 	const double st = sin(at.theta);