Skip to content
GitLab
Commit b1434297 authored by Thomas White's avatar Thomas White
Browse files

Symmetry stuff

parent 012073a3
Hide whitespace changes
Inline Side-by-side
Makefile.am
View file @ b1434297
......@@ -9,11 +9,12 @@ bin_PROGRAMS = src/pattern_sim src/process_hkl src/get_hkl src/indexamajig \
contrib/alter_stream
noinst_PROGRAMS = tests/list_check tests/integration_check \
tests/pr_gradient_check
tests/pr_gradient_check tests/symmetry_check
TESTS = tests/list_check tests/first_merge_check tests/second_merge_check \
tests/third_merge_check tests/fourth_merge_check \
tests/integration_check tests/pr_gradient_check
tests/integration_check tests/pr_gradient_check \
tests/symmetry_check
EXTRA_DIST += tests/first_merge_check tests/second_merge_check \
tests/third_merge_check tests/fourth_merge_check \
......@@ -146,6 +147,10 @@ tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \
src/thread-pool.c src/utils.c src/image.c \
src/hdf5-file.c
tests_symmetry_check_SOURCES = tests/symmetry_check.c src/symmetry.c \
src/utils.c src/thread-pool.c
tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c src/detector.c \
src/cell.c src/geometry.c src/reflist.c \
src/thread-pool.c src/utils.c src/peaks.c \
......
Makefile.in
View file @ b1434297
......@@ -44,11 +44,13 @@ bin_PROGRAMS = src/pattern_sim$(EXEEXT) src/process_hkl$(EXEEXT) \
contrib/alter_stream$(EXEEXT) $(am__EXEEXT_1) $(am__EXEEXT_2)
noinst_PROGRAMS = tests/list_check$(EXEEXT) \
tests/integration_check$(EXEEXT) \
tests/pr_gradient_check$(EXEEXT) $(am__EXEEXT_3)
tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \
$(am__EXEEXT_3)
TESTS = tests/list_check$(EXEEXT) tests/first_merge_check \
tests/second_merge_check tests/third_merge_check \
tests/fourth_merge_check tests/integration_check$(EXEEXT) \
tests/pr_gradient_check$(EXEEXT) $(am__EXEEXT_3)
tests/pr_gradient_check$(EXEEXT) tests/symmetry_check$(EXEEXT) \
$(am__EXEEXT_3)
@BUILD_HDFSEE_TRUE@am__append_1 = src/hdfsee
@BUILD_CUBEIT_TRUE@am__append_2 = src/cubeit
@HAVE_OPENCL_TRUE@am__append_3 = src/diffraction-gpu.c src/cl-utils.c
......@@ -306,6 +308,12 @@ tests_pr_gradient_check_OBJECTS = \
$(am_tests_pr_gradient_check_OBJECTS)
tests_pr_gradient_check_LDADD = $(LDADD)
tests_pr_gradient_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
am_tests_symmetry_check_OBJECTS = tests/symmetry_check.$(OBJEXT) \
src/symmetry.$(OBJEXT) src/utils.$(OBJEXT) \
src/thread-pool.$(OBJEXT)
tests_symmetry_check_OBJECTS = $(am_tests_symmetry_check_OBJECTS)
tests_symmetry_check_LDADD = $(LDADD)
tests_symmetry_check_DEPENDENCIES = $(top_builddir)/lib/libgnu.a
DEFAULT_INCLUDES = -I.@am__isrc@
depcomp = $(SHELL) $(top_srcdir)/build-aux/depcomp
am__depfiles_maybe = depfiles
......@@ -336,7 +344,8 @@ SOURCES = $(contrib_alter_stream_SOURCES) \
$(src_render_hkl_SOURCES) $(src_sum_stack_SOURCES) \
$(tests_gpu_sim_check_SOURCES) \
$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
$(tests_pr_gradient_check_SOURCES)
$(tests_pr_gradient_check_SOURCES) \
$(tests_symmetry_check_SOURCES)
DIST_SOURCES = $(contrib_alter_stream_SOURCES) \
$(src_calibrate_detector_SOURCES) $(src_check_hkl_SOURCES) \
$(src_compare_hkl_SOURCES) $(am__src_cubeit_SOURCES_DIST) \
......@@ -348,7 +357,8 @@ DIST_SOURCES = $(contrib_alter_stream_SOURCES) \
$(src_sum_stack_SOURCES) \
$(am__tests_gpu_sim_check_SOURCES_DIST) \
$(tests_integration_check_SOURCES) $(tests_list_check_SOURCES) \
$(tests_pr_gradient_check_SOURCES)
$(tests_pr_gradient_check_SOURCES) \
$(tests_symmetry_check_SOURCES)
RECURSIVE_TARGETS = all-recursive check-recursive dvi-recursive \
html-recursive info-recursive install-data-recursive \
install-dvi-recursive install-exec-recursive \
......@@ -751,6 +761,9 @@ tests_integration_check_SOURCES = tests/integration_check.c src/peaks.c \
src/thread-pool.c src/utils.c src/image.c \
src/hdf5-file.c
tests_symmetry_check_SOURCES = tests/symmetry_check.c src/symmetry.c \
src/utils.c src/thread-pool.c
tests_pr_gradient_check_SOURCES = tests/pr_gradient_check.c src/detector.c \
src/cell.c src/geometry.c src/reflist.c \
src/thread-pool.c src/utils.c src/peaks.c \
......@@ -1040,6 +1053,11 @@ tests/pr_gradient_check.$(OBJEXT): tests/$(am__dirstamp) \
tests/pr_gradient_check$(EXEEXT): $(tests_pr_gradient_check_OBJECTS) $(tests_pr_gradient_check_DEPENDENCIES) tests/$(am__dirstamp)
@rm -f tests/pr_gradient_check$(EXEEXT)
$(AM_V_CCLD)$(LINK) $(tests_pr_gradient_check_OBJECTS) $(tests_pr_gradient_check_LDADD) $(LIBS)
tests/symmetry_check.$(OBJEXT): tests/$(am__dirstamp) \
tests/$(DEPDIR)/$(am__dirstamp)
tests/symmetry_check$(EXEEXT): $(tests_symmetry_check_OBJECTS) $(tests_symmetry_check_DEPENDENCIES) tests/$(am__dirstamp)
@rm -f tests/symmetry_check$(EXEEXT)
$(AM_V_CCLD)$(LINK) $(tests_symmetry_check_OBJECTS) $(tests_symmetry_check_LDADD) $(LIBS)
mostlyclean-compile:
-rm -f *.$(OBJEXT)
......@@ -1089,6 +1107,7 @@ mostlyclean-compile:
-rm -f tests/integration_check.$(OBJEXT)
-rm -f tests/list_check.$(OBJEXT)
-rm -f tests/pr_gradient_check.$(OBJEXT)
-rm -f tests/symmetry_check.$(OBJEXT)
distclean-compile:
-rm -f *.tab.c
......@@ -1139,6 +1158,7 @@ distclean-compile:
@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/integration_check.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/list_check.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/pr_gradient_check.Po@am__quote@
@AMDEP_TRUE@@am__include@ @am__quote@tests/$(DEPDIR)/symmetry_check.Po@am__quote@
.c.o:
@am__fastdepCC_TRUE@ $(AM_V_CC)depbase=`echo $@ | sed 's|[^/]*$$|$(DEPDIR)/&|;s|\.o$$||'`;\
......
src/diffraction.c
View file @ b1434297
......@@ -100,14 +100,14 @@ static double sym_lookup_intensity(const double *intensities,
int i;
double ret = 0.0;
for ( i=0; i<num_equivs(sym); i++ ) {
for ( i=0; i<num_equivs(sym, NULL); i++ ) {
signed int he;
signed int ke;
signed int le;
double f, val;
get_equiv(sym, i, h, k, l, &he, &ke, &le);
get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
f = (double)lookup_flag(flags, he, ke, le);
val = lookup_intensity(intensities, he, ke, le);
......@@ -127,14 +127,14 @@ static double sym_lookup_phase(const double *phases,
int i;
double ret = 0.0;
for ( i=0; i<num_equivs(sym); i++ ) {
for ( i=0; i<num_equivs(sym, NULL); i++ ) {
signed int he;
signed int ke;
signed int le;
double f, val;
get_equiv(sym, i, h, k, l, &he, &ke, &le);
get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
f = (double)lookup_flag(flags, he, ke, le);
val = lookup_phase(phases, he, ke, le);
......
src/get_hkl.c
View file @ b1434297
......@@ -146,15 +146,17 @@ static RefList *twin_reflections(RefList *in,
Reflection *refl;
RefListIterator *iter;
RefList *out;
SymOpMask *m;
/* FIXME: Check properly by coset decomposition */
if ( num_equivs(holo) < num_equivs(mero) ) {
if ( num_equivs(holo, NULL) < num_equivs(mero, NULL) ) {
ERROR("%s is not a subgroup of %s!\n", symmetry_name(mero),
symmetry_name(holo));
return NULL;
}
out = reflist_new();
m = new_symopmask(holo);
for ( refl = first_refl(in, &iter);
refl != NULL;
......@@ -178,13 +180,15 @@ static RefList *twin_reflections(RefList *in,
total = 0.0;
sigma = 0.0;
skip = 0;
n = num_equivs(holo);
special_position(holo, m, h, k, l);
n = num_equivs(holo, m);
for ( j=0; j<n; j++ ) {
signed int he, ke, le;
signed int hu, ku, lu;
get_equiv(holo, j, h, k, l, &he, &ke, &le);
get_equiv(holo, m, j, h, k, l, &he, &ke, &le);
/* Do we have this reflection?
* We might not have the particular (merohedral)
......@@ -229,15 +233,17 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
Reflection *refl;
RefListIterator *iter;
RefList *out;
SymOpMask *m;
/* FIXME: Check properly */
if ( num_equivs(target) > num_equivs(initial) ) {
if ( num_equivs(target, NULL) > num_equivs(initial, NULL) ) {
ERROR("%s is not a subgroup of %s!\n", symmetry_name(initial),
symmetry_name(target));
return NULL;
}
out = reflist_new();
m = new_symopmask(initial);
for ( refl = first_refl(in, &iter);
refl != NULL;
......@@ -250,7 +256,8 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
get_indices(refl, &h, &k, &l);
intensity = get_intensity(refl);
n = num_equivs(initial);
special_position(initial, m, h, k, l);
n = num_equivs(initial, m);
/* For each equivalent in the higher symmetry group */
for ( j=0; j<n; j++ ) {
......@@ -259,7 +266,7 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
Reflection *new;
/* Get the equivalent */
get_equiv(initial, j, h, k, l, &he, &ke, &le);
get_equiv(initial, m, j, h, k, l, &he, &ke, &le);
/* Put it into the asymmetric unit for the target */
get_asymm(target, he, ke, le, &he, &ke, &le);
......@@ -272,6 +279,8 @@ static RefList *expand_reflections(RefList *in, const SymOpList *target,
}
free_symopmask(m);
return out;
}
......@@ -460,6 +469,9 @@ int main(int argc, char *argv[])
Reflection *refl;
RefListIterator *iter;
SymOpMask *m;
m = new_symopmask(mero);
for ( refl = first_refl(input, &iter);
refl != NULL;
......@@ -471,10 +483,13 @@ int main(int argc, char *argv[])
get_indices(refl, &h, &k, &l);
inty = get_intensity(refl);
inty *= (double)num_equivs(mero);
special_position(mero, m, h, k, l);
inty *= (double)num_equivs(mero, m);
set_int(refl, inty);
}
free_symopmask(m);
}
if ( template ) {
......
src/povray.c
View file @ b1434297
......@@ -42,6 +42,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
int i;
Reflection *refl;
RefListIterator *iter;
SymOpMask *m;
if ( (nproc > MAX_PROC) || (nproc < 1) ) {
ERROR("Number of processes must be a number between 1 and %i\n",
......@@ -167,6 +168,8 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
fprintf(fh, "}\n");
m = new_symopmask(sym);
max = 0.0;
for ( refl = first_refl(list, &iter);
refl != NULL;
......@@ -174,10 +177,9 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
float val;
signed int h, k, l;
SymOpList *sp;
get_indices(refl, &h, &k, &l);
sp = special_position(sym, h, k, l);
special_position(sym, m, h, k, l);
switch ( wght ) {
case WGHT_I :
......@@ -189,7 +191,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
break;
case WGHT_COUNTS :
val = get_redundancy(refl);
val /= (double)num_equivs(sp);
val /= (double)num_equivs(sym, m);
break;
case WGHT_RAWCOUNTS :
val = get_redundancy(refl);
......@@ -201,8 +203,6 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
if ( val > max ) max = val;
free_symoplist(sp);
}
max /= boost;
......@@ -211,7 +211,6 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
max = scale_top;
}
for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
......@@ -220,12 +219,11 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
int s;
float val, p, r, g, b, trans;
int j;
SymOpList *sp;
int neq;
get_indices(refl, &h, &k, &l);
sp = special_position(sym, h, k, l);
neq = num_equivs(sp);
special_position(sym, m, h, k, l);
neq = num_equivs(sym, m);
switch ( wght ) {
case WGHT_I :
......@@ -294,7 +292,7 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
signed int he, ke, le;
float x, y, z;
get_equiv(sp, j, h, k, l, &he, &ke, &le);
get_equiv(sym, m, j, h, k, l, &he, &ke, &le);
x = asx*he + bsx*ke + csx*le;
y = asy*he + bsy*ke + csy*le;
......@@ -311,10 +309,10 @@ int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc,
}
free_symoplist(sp);
}
free_symopmask(m);
fprintf(fh, "\n");
fclose(fh);
......
src/powder_plot.c
View file @ b1434297
......@@ -318,6 +318,9 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
int h, k, l, redundancy;
double q, intensity;
unsigned int nref;
SymOpMask *m;
m = new_symopmask(sym);
nref = num_reflections(image->reflections);
......@@ -330,9 +333,8 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
if ( use_redundancy ) {
redundancy = get_redundancy(refl);
} else {
SymOpList *sp = special_position(sym, h, k, l);
redundancy = num_equivs(sp);
free_symoplist(sp);
special_position(sym, m, h, k, l);
redundancy = num_equivs(sym, m);
}
/* Multiply by 2 to get 1/d (in m^-1) */
......@@ -348,6 +350,8 @@ static unsigned int process_hkl(struct image *image, const SymOpList *sym,
}
free_symopmask(m);
return n_peaks;
}
......
src/process_hkl.c
View file @ b1434297
......@@ -521,6 +521,7 @@ int main(int argc, char *argv[])
if ( sym_str == NULL ) sym_str = strdup("1");
sym = get_pointgroup(sym_str);
STATUS("%s -> %p\n", sym_str, sym);
free(sym_str);
/* Open the data stream */
......@@ -553,7 +554,7 @@ int main(int argc, char *argv[])
ERROR("Invalid indices for '--histogram'\n");
return 1;
}
space_for_hist = n_total_patterns * num_equivs(sym);
space_for_hist = n_total_patterns * num_equivs(sym, NULL);
hist_vals = malloc(space_for_hist * sizeof(double));
free(histo);
STATUS("Histogramming %i %i %i -> ", hist_h, hist_k, hist_l);
......
src/reflist-utils.c
View file @ b1434297
......@@ -109,8 +109,14 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
unsigned char *flags;
Reflection *refl;
RefListIterator *iter;
SymOpMask *mask;
flags = flags_from_list(list);
mask = new_symopmask(sym);
if ( mask == NULL ) {
ERROR("Couldn't create mask for list symmetry check.\n");
return 1;
}
for ( refl = first_refl(list, &iter);
refl != NULL;
......@@ -119,13 +125,19 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
int j;
int found = 0;
signed int h, k, l;
int n;
get_indices(refl, &h, &k, &l);
for ( j=0; j<num_equivs(sym); j++ ) {
special_position(sym, mask, h, k, l);
n = num_equivs(sym, mask);
STATUS("%i equivs: %3i %3i %3i\n", n, h, k, l);
for ( j=0; j<n; j++ ) {
signed int he, ke, le;
get_equiv(sym, j, h, k, l, &he, &ke, &le);
get_equiv(sym, mask, j, h, k, l, &he, &ke, &le);
STATUS("%3i: %3i %3i %3i\n", j, he, ke, le);
if ( abs(he) > INDMAX ) continue;
if ( abs(le) > INDMAX ) continue;
......@@ -137,12 +149,15 @@ int check_list_symmetry(RefList *list, const SymOpList *sym)
if ( found > 1 ) {
free(flags);
free_symopmask(mask);
STATUS("Found %i %i %i twice\n", h, k, l);
return 1; /* Symmetry is wrong! */
}
}
free(flags);
free_symopmask(mask);
return 0;
}
......@@ -155,11 +170,11 @@ int find_equiv_in_list(RefList *list, signed int h, signed int k,
int i;
int found = 0;
for ( i=0; i<num_equivs( sym); i++ ) {
for ( i=0; i<num_equivs(sym, NULL); i++ ) {
signed int he, ke, le;
Reflection *f;
get_equiv(sym, i, h, k, l, &he, &ke, &le);
get_equiv(sym, NULL, i, h, k, l, &he, &ke, &le);
f = find_refl(list, he, ke, le);
/* There must only be one equivalent. If there are more, it
......
src/render_hkl.c
View file @ b1434297
......@@ -96,12 +96,15 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
{
Reflection *refl;
RefListIterator *iter;
SymOpMask *m;
if ( dctx == NULL ) {
*max_u = 0.0; *max_v = 0.0; *max_val = 0.0;
*max_res = 0.0; *max_ux = 0; *max_uy = 0;
}
m = new_symopmask(sym);
/* Iterate over all reflections */
for ( refl = first_refl(list, &iter);
refl != NULL;
......@@ -111,18 +114,17 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
signed int ha, ka, la;
int xi, yi;
int i, n;
SymOpList *sp;
get_indices(refl, &ha, &ka, &la);
sp = special_position(sym, ha, ka, la);
n = num_equivs(sp);
special_position(sym, m, ha, ka, la);
n = num_equivs(sym, m);
for ( i=0; i<n; i++ ) {
signed int h, k, l;
get_equiv(sp, i, ha, ka, la, &h, &k, &l);
get_equiv(sym, m, i, ha, ka, la, &h, &k, &l);
/* Is the reflection in the zone? */
if ( h*zh + k*zk + l*zl != 0 ) continue;
......@@ -194,9 +196,9 @@ static void draw_circles(signed int xh, signed int xk, signed int xl,
}
free_symoplist(sp);
}
free_symopmask(m);
}
......
src/symmetry.c
View file @ b1434297
......@@ -37,19 +37,6 @@
*/
enum lattice_type
{
L_TRICLINIC,
L_MONOCLINIC,
L_ORTHORHOMBIC,
L_TETRAGONAL,
L_RHOMBOHEDRAL,
L_TRIGONAL,
L_HEXAGONAL,
L_CUBIC,
};
struct sym_op
{
signed int *h;
......@@ -84,6 +71,13 @@ struct _symoplist
};
struct _symopmask
{
const SymOpList *list;
int *mask;
};
static void alloc_ops(SymOpList *ops)
{
......@@ -92,6 +86,35 @@ static void alloc_ops(SymOpList *ops)
}
/**
* new_symopmask:
*
* Returns: a new %SymOpMask, which you can use when filtering out special
* reflections.
**/
SymOpMask *new_symopmask(const SymOpList *list)
{
SymOpMask *m;
int i;
m = malloc(sizeof(struct _symopmask));
if ( m == NULL ) return NULL;
m->list = list;
m->mask = malloc(sizeof(int)*list->n_ops);
if ( m->mask == NULL ) {
free(m);
return NULL;
}
for ( i=0; i<list->n_ops; i++ ) {
m->mask[i] = 1;
}
return m;
}
/* Creates a new SymOpList */
static SymOpList *new_symoplist()
{
......@@ -101,6 +124,7 @@ static SymOpList *new_symoplist()
new->max_ops = 16;
new->n_ops = 0;
new->ops = NULL;
new->divisors = NULL;
new->name = NULL;
new->num_equivs = 1;
alloc_ops(new);
......@@ -128,88 +152,33 @@ void free_symoplist(SymOpList *ops)
free(ops);
}
static int is_identity(signed int *h, signed int *k, signed int *l)
/**
* free_symopmask:
*
* Frees a %SymOpMask and all associated resources.
**/
void free_symopmask(SymOpMask *m)
{
if ( (h[0]!=1) || (h[1]!=0) || (h[2]!=0) ) return 0;
if ( (k[0]!=0) || (k[1]!=1) || (k[2]!=0) ) return 0;
if ( (l[0]!=0) || (l[1]!=0) || (l[2]!=1) ) return 0;
return 1;
if ( m == NULL ) return;
free(m->mask);
free(m);
}
/* Calculate the order of the operation "M", which is the lowest
* integer n such that M^n = I. */
static int order_of_op(signed int *hin, signed int *kin, signed int *lin)
{
int n;
signed int h[3];
signed int k[3];
signed int l[3];
memcpy(h, hin, 3*sizeof(signed int));
memcpy(k, kin, 3*sizeof(signed int));
memcpy(l, lin, 3*sizeof(signed int));
for ( n=1; n<6; n++ ) {
signed int hnew[3];
signed int knew[3];
signed int lnew[3];
/* Yay matrices */
hnew[0] = h[0]*h[0] + h[1]*k[0] + h[2]*l[0];
hnew[1] = h[0]*h[1] + h[1]*k[1] + h[2]*l[1];
hnew[2] = h[0]*h[2] + h[1]*k[2] + h[2]*l[2];
knew[0] = k[0]*h[0] + k[1]*k[0] + k[2]*l[0];
knew[1] = k[0]*h[1] + k[1]*k[1] + k[2]*l[1];
knew[2] = k[0]*h[2] + k[1]*k[2] + k[2]*l[2];
lnew[0] = l[0]*h[0] + l[1]*