diff options
Diffstat (limited to 'libcrystfel')
| -rw-r--r-- | libcrystfel/Makefile.am | 5 | ||||
| -rw-r--r-- | libcrystfel/src/integer_matrix.c | 386 | ||||
| -rw-r--r-- | libcrystfel/src/integer_matrix.h | 68 | ||||
| -rw-r--r-- | libcrystfel/src/symmetry.c | 265 |
4 files changed, 684 insertions, 40 deletions
diff --git a/libcrystfel/Makefile.am b/libcrystfel/Makefile.am index 28d483ac..2a02a20b 100644 --- a/libcrystfel/Makefile.am +++ b/libcrystfel/Makefile.am @@ -8,7 +8,7 @@ libcrystfel_la_SOURCES = src/reflist.c src/utils.c src/cell.c src/detector.c \ src/symmetry.c src/stream.c src/peaks.c \ src/reflist-utils.c src/filters.c \ src/render.c src/index.c src/dirax.c src/mosflm.c \ - src/cell-utils.c + src/cell-utils.c src/integer_matrix.c if HAVE_FFTW libcrystfel_la_SOURCES += src/reax.c @@ -23,7 +23,8 @@ libcrystfel_la_include_HEADERS = src/beam-parameters.h src/hdf5-file.h \ src/geometry.h src/peaks.h src/stream.h \ src/render.h src/index.h src/image.h \ src/filters.h src/dirax.h src/mosflm.h \ - src/index-priv.h src/reax.h src/cell-utils.h + src/index-priv.h src/reax.h src/cell-utils.h \ + src/integer_matrix.h INCLUDES = "-I$(top_srcdir)/data" AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -Wall diff --git a/libcrystfel/src/integer_matrix.c b/libcrystfel/src/integer_matrix.c new file mode 100644 index 00000000..5aefb958 --- /dev/null +++ b/libcrystfel/src/integer_matrix.c @@ -0,0 +1,386 @@ +/* + * integer_matrix.c + * + * A small integer matrix library + * + * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +#include <stdlib.h> +#include <stdio.h> +#include <string.h> +#include <assert.h> + +#include "integer_matrix.h" + + +/** + * SECTION:integer_matrix + * @short_description: Integer matrices + * @title: Integer matrices + * @section_id: + * @see_also: + * @include: "integer_matrix.h" + * @Image: + * + * Routines to handle integer matrix + */ + + +struct _integermatrix +{ + unsigned int rows; + unsigned int cols; + + signed int *v; +}; + + +/** + * intmat_new: + * + * Allocates a new %IntegerMatrix with all elements set to zero. + * + * Returns: a new %IntegerMatrix, or NULL on error. + **/ +IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols) +{ + IntegerMatrix *m; + + m = malloc(sizeof(IntegerMatrix)); + if ( m == NULL ) return NULL; + + m->v = calloc(rows*cols, sizeof(signed int)); + if ( m->v == NULL ) { + free(m); + return NULL; + } + + m->rows = rows; + m->cols = cols; + + return m; +} + + +/** + * intmat_free: + * @m: An %IntegerMatrix + * + * Frees @m. + **/ +void intmat_free(IntegerMatrix *m) +{ + free(m->v); + free(m); +} + + +/** + * intmat_set: + * @m: An %IntegerMatrix + * @i: row number to set + * @j: column number to set + * @v: value to set to + * + * Sets the @i,@j element of @m to @v. + **/ +void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j, signed int v) +{ + assert(i < m->rows); + assert(j < m->cols); + m->v[j + m->cols*i] = v; +} + + +/** + * intmat_get: + * @m: An %IntegerMatrix + * @i: column number to set + * @j: row number to set + * + * Gets the @i,@j element of @m. + * + * Returns: the @i,@j element of @m. + **/ +signed int intmat_get(IntegerMatrix *m, unsigned int i, unsigned int j) +{ + assert(i < m->rows); + assert(j < m->cols); + return m->v[j + m->cols*i]; +} + + +/** + * intmat_intvec_mult: + * @m: An %IntegerMatrix + * @vec: An array of signed integers + * + * Multiplies the matrix @m by the vector @vec. The size of @vec must equal the + * number of columns in @m, and the size of the result equals the number of rows + * in @m. + * + * Returns: a newly allocated array of signed integers containing the answer, + * or NULL on error. + **/ +signed int *intmat_intvec_mult(IntegerMatrix *m, signed int *vec) +{ + signed int *ans; + unsigned int i; + + ans = malloc(m->rows * sizeof(signed int)); + if ( ans == NULL ) return NULL; + + for ( i=0; i<m->rows; i++ ) { + + unsigned int j; + + ans[i] = 0; + for ( j=0; j<m->cols; j++ ) { + ans[i] += intmat_get(m, i, j) * vec[j]; + } + + } + + return ans; +} + + +/** + * intmat_intmat_mult: + * @a: An %IntegerMatrix + * @b: An %IntegerMatrix + * + * Multiplies the matrix @a by the matrix @b. + * + * Returns: a newly allocated %IntegerMatrix containing the answer, or NULL on + * error. + **/ +IntegerMatrix *intmat_intmat_mult(IntegerMatrix *a, IntegerMatrix *b) +{ + unsigned int i, j; + IntegerMatrix *ans; + + if ( a->cols != b->rows ) return NULL; + + ans = intmat_new(a->rows, a->cols); + if ( ans == NULL ) return NULL; + + for ( i=0; i<ans->rows; i++ ) { + for ( j=0; j<ans->cols; j++ ) { + + unsigned int k; + signed int r = 0; + + for ( k=0; k<a->cols; k++ ) { /* a->cols == b->rows */ + r += intmat_get(a, i, k) * intmat_get(b, k, j); + } + intmat_set(ans, i, j, r); + + } + } + + return ans; +} + + +static IntegerMatrix *delete_row_and_column(IntegerMatrix *m, + unsigned int di, unsigned int dj) +{ + IntegerMatrix *n; + unsigned int i, j; + + n = intmat_new(m->rows-1, m->cols-1); + if ( n == NULL ) return NULL; + + for ( i=0; i<n->rows; i++ ) { + for ( j=0; j<n->cols; j++ ) { + + signed int val; + unsigned int gi, gj; + + gi = (i>=di) ? i+1 : i; + gj = (j>=dj) ? j+1 : j; + val = intmat_get(m, gi, gj); + intmat_set(n, i, j, val); + + } + } + + return n; +} + + +/** + * intmat_det: + * @m: An %IntegerMatrix + * + * Calculates the determinant of @m. Inefficiently. + * + * Returns: the determinant of @m. + **/ +signed int intmat_det(IntegerMatrix *m) +{ + unsigned int i, j; + signed int det = 0; + + assert(m->rows == m->cols); /* Otherwise determinant doesn't exist */ + + if ( m->rows == 2 ) { + return intmat_get(m, 0, 0)*intmat_get(m, 1, 1) + - intmat_get(m, 0, 1)*intmat_get(m, 1, 0); + } + + i = 0; /* Fixed */ + for ( j=0; j<m->cols; j++ ) { + + signed int C; + IntegerMatrix *n; + signed int t; + + n = delete_row_and_column(m, i, j); + if ( n == NULL ) { + fprintf(stderr, "Failed to allocate matrix.\n"); + return 0; + } + + /* -1 if odd, +1 if even */ + t = (i+j) & 0x1 ? -1 : +1; + + C = t * intmat_det(n); + intmat_free(n); + + det += intmat_get(m, i, j) * C; + + } + + return det; +} + + +static IntegerMatrix *intmat_cofactors(IntegerMatrix *m) +{ + IntegerMatrix *n; + signed int i, j; + + n = intmat_new(m->rows, m->cols); + if ( n == NULL ) return NULL; + + for ( i=0; i<n->rows; i++ ) { + for ( j=0; j<n->cols; j++ ) { + + signed int C; + IntegerMatrix *M; + signed int t; + + M = delete_row_and_column(m, i, j); + if ( M == NULL ) { + fprintf(stderr, "Failed to allocate matrix.\n"); + return 0; + } + + /* -1 if odd, +1 if even */ + t = (i+j) & 0x1 ? -1 : +1; + + C = t * intmat_det(M); + intmat_free(M); + + intmat_set(n, i, j, C); + + } + } + + return n; +} + + +/** + * intmat_inverse: + * @m: An %IntegerMatrix + * + * Calculates the inverse of @m. Inefficiently. + * + * Returns: the inverse of @m, or NULL on error. + **/ +IntegerMatrix *intmat_inverse(IntegerMatrix *m) +{ + IntegerMatrix *adjugateT; + IntegerMatrix *inverse; + unsigned int i, j; + signed int det; + + det = intmat_det(m); + if ( (det != +1) && (det != -1) ) { + fprintf(stderr, "Inverse matrix not an integer matrix.\n"); + return NULL; + } + + adjugateT = intmat_cofactors(m); + if ( adjugateT == NULL ) return NULL; + + inverse = intmat_new(m->cols, m->rows); /* The other way round */ + if ( inverse == NULL ) return NULL; + + for ( i=0; i<inverse->rows; i++ ) { + for ( j=0; j<inverse->cols; j++ ) { + + signed int v; + + v = intmat_get(adjugateT, j, i); + + /* 1/-1 = -1 and 1/+1 = +1, and these are the only two cases */ + intmat_set(inverse, i, j, v*det); + + } + } + + intmat_free(adjugateT); + + return inverse; +} + + +/** + * intmat_print + * @m: An %IntegerMatrix + * + * Prints @m to stderr. + * + */ +void intmat_print(IntegerMatrix *m) +{ + unsigned int i, j; + + for ( i=0; i<m->rows; i++ ) { + + fprintf(stderr, "[ "); + for ( j=0; j<m->cols; j++ ) { + fprintf(stderr, "%4i ", intmat_get(m, i, j)); + } + fprintf(stderr, "]\n"); + } +} diff --git a/libcrystfel/src/integer_matrix.h b/libcrystfel/src/integer_matrix.h new file mode 100644 index 00000000..41780629 --- /dev/null +++ b/libcrystfel/src/integer_matrix.h @@ -0,0 +1,68 @@ +/* + * integer_matrix.h + * + * A small integer matrix library + * + * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. + * + * Authors: + * 2012 Thomas White <taw@physics.org> + * + * This file is part of CrystFEL. + * + * CrystFEL is free software: you can redistribute it and/or modify + * it under the terms of the GNU General Public License as published by + * the Free Software Foundation, either version 3 of the License, or + * (at your option) any later version. + * + * CrystFEL is distributed in the hope that it will be useful, + * but WITHOUT ANY WARRANTY; without even the implied warranty of + * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the + * GNU General Public License for more details. + * + * You should have received a copy of the GNU General Public License + * along with CrystFEL. If not, see <http://www.gnu.org/licenses/>. + * + */ + +#ifndef INTEGER_MATRIX_H +#define INTEGER_MATRIX_H + +#ifdef HAVE_CONFIG_H +#include <config.h> +#endif + +/** + * IntegerMatrix + * + * The IntegerMatrix is an opaque data structure representing an integer matrix. + **/ +typedef struct _integermatrix IntegerMatrix; + + +/* Alloc/dealloc */ +extern IntegerMatrix *intmat_new(unsigned int rows, unsigned int cols); +extern void intmat_free(IntegerMatrix *m); + +/* Get/set */ +extern void intmat_set(IntegerMatrix *m, unsigned int i, unsigned int j, + signed int v); +extern signed int intmat_get(IntegerMatrix *m, unsigned int i, unsigned int j); + +/* Matrix-(int)vector multiplication */ +extern signed int *intmat_intvec_mult(IntegerMatrix *m, signed int *vec); + +/* Matrix-matrix multiplication */ +extern IntegerMatrix *intmat_intmat_mult(IntegerMatrix *a, IntegerMatrix *b); + +/* Inverse */ +extern IntegerMatrix *intmat_inverse(IntegerMatrix *m); + +/* Determinant */ +extern signed int intmat_det(IntegerMatrix *m); + +/* Diagnostics */ +extern void intmat_print(IntegerMatrix *m); + +#endif /* INTEGER_MATRIX_H */ diff --git a/libcrystfel/src/symmetry.c b/libcrystfel/src/symmetry.c index 52a9aae6..14681434 100644 --- a/libcrystfel/src/symmetry.c +++ b/libcrystfel/src/symmetry.c @@ -38,6 +38,7 @@ #include "symmetry.h" #include "utils.h" +#include "integer_matrix.h" /** @@ -57,7 +58,7 @@ struct sym_op { signed int *h; signed int *k; - signed int *l; /* Contributions to h, k and l from h, k, i and l */ + signed int *l; /* Contributions to h, k and l from h, k and l */ int order; }; @@ -359,6 +360,100 @@ static SymOpList *expand_ops(SymOpList *s) } +/* Transform all the operations in a SymOpList by a given matrix. + * The matrix must have a determinant of +/- 1 (otherwise its inverse would + * not also be an integer matrix). */ +static void transform_ops(SymOpList *s, signed int *na, + signed int *nb, + signed int *nc) +{ + int n, i; + IntegerMatrix *t, *inv; + signed int det; + + t = intmat_new(3, 3); + if ( t == NULL ) { + ERROR("Failed to allocate matrix.\n"); + return; + } + + intmat_set(t, 0, 0, na[0]); + intmat_set(t, 1, 0, na[1]); + intmat_set(t, 2, 0, na[2]); + intmat_set(t, 0, 1, nb[0]); + intmat_set(t, 1, 1, nb[1]); + intmat_set(t, 2, 1, nb[2]); + intmat_set(t, 0, 2, nc[0]); + intmat_set(t, 1, 2, nc[1]); + intmat_set(t, 2, 2, nc[2]); + + det = intmat_det(t); + if ( det == -1 ) { + ERROR("Warning: mirrored SymOpList.\n"); + } else if ( det != 1 ) { + ERROR("Invalid transformation for SymOpList.\n"); + return; + } + + inv = intmat_inverse(t); + if ( inv == NULL ) { + ERROR("Failed to invert matrix.\n"); + return; + } + + n = num_ops(s); + for ( i=0; i<n; i++ ) { + + IntegerMatrix *m, *r, *f; + + m = intmat_new(3, 3); + if ( m == NULL ) { + ERROR("Failed to allocate matrix.\n"); + return; + } + + intmat_set(m, 0, 0, s->ops[i].h[0]); + intmat_set(m, 1, 0, s->ops[i].h[1]); + intmat_set(m, 2, 0, s->ops[i].h[2]); + intmat_set(m, 0, 1, s->ops[i].k[0]); + intmat_set(m, 1, 1, s->ops[i].k[1]); + intmat_set(m, 2, 1, s->ops[i].k[2]); + intmat_set(m, 0, 2, s->ops[i].l[0]); + intmat_set(m, 1, 2, s->ops[i].l[1]); + intmat_set(m, 2, 2, s->ops[i].l[2]); + + r = intmat_intmat_mult(m, t); + if ( r == NULL ) { + ERROR("Matrix multiplication failed.\n"); + return; + } + intmat_free(m); + + f = intmat_intmat_mult(inv, r); + if ( f == NULL ) { + ERROR("Matrix multiplication failed.\n"); + return; + } + intmat_free(r); + + s->ops[i].h[0] = intmat_get(f, 0, 0); + s->ops[i].h[1] = intmat_get(f, 1, 0); + s->ops[i].h[2] = intmat_get(f, 2, 0); + s->ops[i].k[0] = intmat_get(f, 0, 1); + s->ops[i].k[1] = intmat_get(f, 1, 1); + s->ops[i].k[2] = intmat_get(f, 2, 1); + s->ops[i].l[0] = intmat_get(f, 0, 2); + s->ops[i].l[1] = intmat_get(f, 1, 2); + s->ops[i].l[2] = intmat_get(f, 2, 2); + intmat_free(f); + + } + + intmat_free(t); + intmat_free(inv); +} + + /********************************* Triclinic **********************************/ static SymOpList *make_1bar() @@ -390,15 +485,6 @@ static SymOpList *make_2m() } -static SymOpList *make_2_uab() -{ - SymOpList *new = new_symoplist(); - add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,-1), 2); /* 2 // k */ - new->name = strdup("2_uab"); - return expand_ops(new); -} - - static SymOpList *make_2() { SymOpList *new = new_symoplist(); @@ -531,17 +617,6 @@ static SymOpList *make_4mmm() } -static SymOpList *make_4mmm_uaa() -{ - SymOpList *new = new_symoplist(); - add_symop(new, v(1,0,0,0), v(0,0,0,1), v(0,-1,0,0), 4); /* 4 // h */ - add_symop(new, v(1,0,0,0), v(0,-1,0,0), v(0,0,0,1), 2); /* m -| k */ - add_symop(new, v(-1,0,0,0), v(0,1,0,0), v(0,0,0,1), 2); /* m -| h */ - new->name = strdup("4/mmm_uaa"); - return expand_ops(new); -} - - /************************** Trigonal (Rhombohedral) ***************************/ static SymOpList *make_3_R() @@ -817,21 +892,7 @@ static SymOpList *make_m3barm() } -/** - * get_pointgroup: - * @sym: A string representation of a point group - * - * This function parses @sym and returns the corresponding %SymOpList. - * In the string representation of the point group, use a preceding minus sign - * for any character which would have a "bar". Trigonal groups must be suffixed - * with either "_H" or "_R" for a hexagonal or rhombohedral lattice - * respectively. - * - * Examples: -1 1 2/m 2 m mmm 222 mm2 4/m 4 -4 4/mmm 422 -42m -4m2 4mm - * 3_R -3_R 32_R 3m_R -3m_R 3_H -3_H 321_H 312_H 3m1_H 31m_H -3m1_H -31m_H - * 6/m 6 -6 6/mmm 622 -62m -6m2 6mm 23 m-3 432 -43m m-3m. - **/ -SymOpList *get_pointgroup(const char *sym) +static SymOpList *getpg_old(const char *sym) { /* Triclinic */ if ( strcmp(sym, "-1") == 0 ) return make_1bar(); @@ -840,7 +901,6 @@ SymOpList *get_pointgroup(const char *sym) /* Monoclinic */ if ( strcmp(sym, "2/m") == 0 ) return make_2m(); if ( strcmp(sym, "2") == 0 ) return make_2(); - if ( strcmp(sym, "2_uab") == 0 ) return make_2_uab(); if ( strcmp(sym, "m") == 0 ) return make_m(); /* Orthorhombic */ @@ -853,7 +913,6 @@ SymOpList *get_pointgroup(const char *sym) if ( strcmp(sym, "4") == 0 ) return make_4(); if ( strcmp(sym, "-4") == 0 ) return make_4bar(); if ( strcmp(sym, "4/mmm") == 0 ) return make_4mmm(); - if ( strcmp(sym, "4/mmm_uaa") == 0 ) return make_4mmm_uaa(); if ( strcmp(sym, "422") == 0 ) return make_422(); if ( strcmp(sym, "-42m") == 0 ) return make_4bar2m(); if ( strcmp(sym, "-4m2") == 0 ) return make_4barm2(); @@ -898,6 +957,136 @@ SymOpList *get_pointgroup(const char *sym) } +static int char_count(const char *a, char b) +{ + size_t i; + int n; + + i = 0; n = 0; + do { + if ( a[i] == b ) n++; + if ( a[i] == '\0' ) return n; + i++; + } while ( 1 ); +} + + +static SymOpList *getpg_old_ua(const char *sym, size_t s) +{ + char ua; + char *pg_type; + SymOpList *pg; + + if ( strncmp(sym+s, "ua", 2) == 0 ) { + ua = sym[s+2]; + } else { + ERROR("Unrecognised point group '%s'\n", sym); + return NULL; + } + + pg_type = strndup(sym, s-1); + if ( pg_type == NULL ) { + ERROR("Couldn't allocate string.\n"); + return NULL; + } + + pg = getpg_old(pg_type); + if ( pg == NULL ) { + ERROR("Unrecognised point group type '%s'\n", + pg_type); + return NULL; + } + free(pg_type); + + switch ( ua ) { + + case 'a' : + transform_ops(pg, v(0,0,0,1), + v(0,1,0,0), + v(-1,0,0,0)); + break; + + case 'b' : + transform_ops(pg, v(1,0,0,0), + v(0,0,0,1), + v(0,-1,0,0)); + break; + + case 'c' : + /* No transformation needed */ + break; + + default : + ERROR("Bad unique axis '%c'\n", ua); + free_symoplist(pg); + return NULL; + + } + + return pg; +} + + +/** + * get_pointgroup: + * @sym: A string representation of a point group + * + * This function parses @sym and returns the corresponding %SymOpList. + * In the string representation of the point group, use a preceding minus sign + * for any character which would have a "bar". Trigonal groups must be suffixed + * with either "_H" or "_R" for a hexagonal or rhombohedral lattice + * respectively. + * + * Examples: -1 1 2/m 2 m mmm 222 mm2 4/m 4 -4 4/mmm 422 -42m -4m2 4mm + * 3_R -3_R 32_R 3m_R -3m_R 3_H -3_H 321_H 312_H 3m1_H 31m_H -3m1_H -31m_H + * 6/m 6 -6 6/mmm 622 -62m -6m2 6mm 23 m-3 432 -43m m-3m. + **/ +SymOpList *get_pointgroup(const char *sym) +{ + int n_space, n_underscore; + + n_space = char_count(sym, ' '); + n_underscore = char_count(sym, '_'); + + if ( n_space == 0 ) { + + /* No spaces nor underscores -> old system */ + if ( n_underscore == 0 ) return getpg_old(sym); + + /* No spaces and 1 underscore -> old system + lattice or UA */ + if ( n_underscore == 1 ) { + + const char *s; + + s = strchr(sym, '_'); + assert(s != NULL); + s++; + + /* Old system with H/R lattice? */ + if ( (s[0] == 'H') || (s[0] == 'R') ) { + return getpg_old(sym); + } + + /* Old system with unique axis */ + return getpg_old_ua(sym, s-sym); + + } + + } + + if ( n_space > 4 ) { + ERROR("Unrecognised point group '%s'\n", sym); + return NULL; + } + + /* 1-4 spaces -> new system, possibly with UA and/or lattice */ + + /* FIXME: Implementation */ + + return NULL; +} + + static void do_op(const struct sym_op *op, signed int h, signed int k, signed int l, signed int *he, signed int *ke, signed int *le) |