diff options
| -rw-r--r-- | src/Makefile.am | 2 | ||||
| -rw-r--r-- | src/Makefile.in | 4 | ||||
| -rw-r--r-- | src/compare_hkl.c | 36 | ||||
| -rw-r--r-- | src/get_hkl.c | 65 | ||||
| -rw-r--r-- | src/symmetry.c | 47 | ||||
| -rw-r--r-- | src/symmetry.h | 4 |
6 files changed, 83 insertions, 75 deletions
diff --git a/src/Makefile.am b/src/Makefile.am index b9b70c22..a61ebb6b 100644 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -39,7 +39,7 @@ get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c symmetry.c get_hkl_LDADD = @LIBS@ compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \ - statistics.c + statistics.c symmetry.c compare_hkl_LDADD = @LIBS@ powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \ diff --git a/src/Makefile.in b/src/Makefile.in index 7c8c1028..7cd5226c 100644 --- a/src/Makefile.in +++ b/src/Makefile.in @@ -60,7 +60,7 @@ calibrate_detector_OBJECTS = $(am_calibrate_detector_OBJECTS) calibrate_detector_DEPENDENCIES = am_compare_hkl_OBJECTS = compare_hkl.$(OBJEXT) sfac.$(OBJEXT) \ cell.$(OBJEXT) utils.$(OBJEXT) reflections.$(OBJEXT) \ - statistics.$(OBJEXT) + statistics.$(OBJEXT) symmetry.$(OBJEXT) compare_hkl_OBJECTS = $(am_compare_hkl_OBJECTS) compare_hkl_DEPENDENCIES = am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \ @@ -256,7 +256,7 @@ indexamajig_LDADD = @LIBS@ get_hkl_SOURCES = get_hkl.c sfac.c cell.c utils.c reflections.c symmetry.c get_hkl_LDADD = @LIBS@ compare_hkl_SOURCES = compare_hkl.c sfac.c cell.c utils.c reflections.c \ - statistics.c + statistics.c symmetry.c compare_hkl_LDADD = @LIBS@ powder_plot_SOURCES = powder_plot.c cell.c utils.c image.c hdf5-file.c \ diff --git a/src/compare_hkl.c b/src/compare_hkl.c index 445798d6..27ca438a 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -25,6 +25,7 @@ #include "sfac.h" #include "reflections.h" #include "statistics.h" +#include "symmetry.h" static void show_help(const char *s) @@ -45,6 +46,7 @@ int main(int argc, char *argv[]) int c; double *ref1; double *ref2; + double *ref2_transformed; double *out; UnitCell *cell; char *outfile = NULL; @@ -64,7 +66,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "ho:", longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "ho:y:", longopts, NULL)) != -1) { switch (c) { case 'h' : @@ -110,36 +112,44 @@ int main(int argc, char *argv[]) return 1; } - /* Find common reflections */ - icommon = intersection_items(i1, i2); - ncom = num_items(icommon); /* List for output scale factor map */ out = new_list_intensity(); - for ( i=0; i<ncom; i++ ) { - double i1, i2; + /* Find common reflections (taking symmetry into account) */ + icommon = new_items(); + ref2_transformed = new_list_intensity(); + for ( i=0; i<num_items(i1); i++ ) { + struct refl_item *it; signed int h, k, l; + signed int he, ke, le; + double val1, val2; - it = get_item(icommon, i); + it = get_item(i1, i); h = it->h; k = it->k; l = it->l; - i1 = lookup_intensity(ref1, h, k, l); - i2 = lookup_intensity(ref2, h, k, l); + if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) { + continue; + } - set_intensity(out, h, k, l, i1/i2); + val1 = lookup_intensity(ref1, h, k, l); + val2 = lookup_intensity(ref2, he, ke, le); + set_intensity(ref2_transformed, h, k, l, val2); + set_intensity(out, h, k, l, val1/val2); + add_item(icommon, h, k, l); } + ncom = num_items(icommon); STATUS("%i,%i reflections: %i in common\n", num_items(i1), num_items(i2), ncom); - R2 = stat_r2(ref1, ref2, icommon, &scale); + R2 = stat_r2(ref1, ref2_transformed, icommon, &scale); STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale); - Rmerge = stat_rmerge(ref1, ref2, icommon, &scale); + Rmerge = stat_rmerge(ref1, ref2_transformed, icommon, &scale); STATUS("Rmerge = %5.4f %% (scale=%5.2e)\n", Rmerge*100.0, scale); - pearson = stat_pearson(ref1, ref2, icommon); + pearson = stat_pearson(ref1, ref2_transformed, icommon); STATUS("Pearson r = %5.4f\n", pearson); if ( outfile != NULL ) { diff --git a/src/get_hkl.c b/src/get_hkl.c index 1052d786..14f79266 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -72,42 +72,6 @@ static void noisify_reflections(double *ref) } -static void scold_user_about_symmetry(signed int h, signed int k, signed int l, - signed int he, signed int ke, - signed int le) -{ - ERROR("Merohedrally equivalent reflection (%i %i %i) found for " - "%i %i %i.\n", he, ke, le, h, k, l); - ERROR("This indicates that you lied to me about the symmetry of the " - "input reflections. "); - ERROR("I won't be able to give you a meaningful result in this " - "situation, so I'm going to give up right now. "); - ERROR("Please reconsider your previous processing of the data, and " - "perhaps try again with a lower symmetry for the '-y' option.\n"); -} - - -static int find_unique_equiv(ReflItemList *items, signed int h, signed int k, - signed int l, const char *mero, signed int *hu, - signed int *ku, signed int *lu) -{ - int i; - - for ( i=0; i<num_equivs(h, k, l, mero); i++ ) { - - signed int he, ke, le; - get_equiv(h, k, l, &he, &ke, &le, mero, i); - if ( find_item(items, he, ke, le) ) { - *hu = he; *ku = ke; *lu = le; - return 1; - } - - } - - return 0; -} - - static ReflItemList *twin_reflections(double *ref, ReflItemList *items, const char *holo, const char *mero) { @@ -125,28 +89,14 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items, int skip; it = get_item(items, i); - h = it->h; k = it->k; l = it->l; - /* None of the equivalent reflections should exist in the - * input dataset. That would indicate that the user lied about - * the input symmetry. - * - * Start from j=1 to ignore the reflection itself. - */ - for ( j=1; j<num_equivs(h, k, l, mero); j++ ) { - - signed int he, ke, le; - get_equiv(h, k, l, &he, &ke, &le, mero, j); - if ( !find_item(items, he, ke, le) ) continue; - - scold_user_about_symmetry(h, k, l, he, ke, le); - abort(); - - } - /* It doesn't matter if the reflection wasn't actually the one - * we define as being in the asymmetric unit cell, as long as - * things aren't confused by there being more than one of it. + /* There is a many-to-one correspondence between reflections + * in the merohedral and holohedral unit cells. Do the + * calculation only once for each reflection in the holohedral + * cell, which contains fewer reflections. */ + get_asymm(it->h, it->k, it->l, &h, &k, &l, holo); + if ( find_item(new, h, k, l) ) continue; n = num_equivs(h, k, l, holo); @@ -163,9 +113,6 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items, * We might not have the particular (merohedral) * equivalent which belongs to our definition of the * asymmetric unit cell, so check them all. - * - * We checked earlier that there's only one of these - * for each reflection. */ if ( !find_unique_equiv(items, he, ke, le, mero, &hu, &ku, &lu) ) { diff --git a/src/symmetry.c b/src/symmetry.c index 4bc8d115..54b7f5cc 100644 --- a/src/symmetry.c +++ b/src/symmetry.c @@ -352,3 +352,50 @@ ReflItemList *get_twins(ReflItemList *items, const char *holo, const char *mero) return ops; } + + +static void scold_user_about_symmetry(signed int h, signed int k, signed int l, + signed int he, signed int ke, + signed int le) +{ + ERROR("Symmetrically equivalent reflection (%i %i %i) found for " + "%i %i %i in the input.\n", he, ke, le, h, k, l); + ERROR("This indicates that you lied to me about the symmetry of the " + "input reflections. "); + ERROR("I won't be able to give you a meaningful result in this " + "situation, so I'm going to give up right now. "); + ERROR("Please reconsider your previous processing of the data, and " + "perhaps try again with a lower symmetry for the '-y' option.\n"); + abort(); +} + + +int find_unique_equiv(ReflItemList *items, signed int h, signed int k, + signed int l, const char *mero, signed int *hu, + signed int *ku, signed int *lu) +{ + int i; + int found = 0; + + for ( i=0; i<num_equivs(h, k, l, mero); i++ ) { + + signed int he, ke, le; + int f; + get_equiv(h, k, l, &he, &ke, &le, mero, i); + f = find_item(items, he, ke, le); + + /* There must only be one equivalent. If there are more, it + * indicates that the user lied about the input symmetry. */ + if ( f && found ) { + scold_user_about_symmetry(he, ke, le, *hu, *ku, *lu); + } + + if ( f && !found ) { + *hu = he; *ku = ke; *lu = le; + found = 1; + } + + } + + return found; +} diff --git a/src/symmetry.h b/src/symmetry.h index 437066ad..bb4ffa60 100644 --- a/src/symmetry.h +++ b/src/symmetry.h @@ -40,5 +40,9 @@ extern const char *get_holohedral(const char *sym); extern ReflItemList *get_twins(ReflItemList *items, const char *holo, const char *mero); +extern int find_unique_equiv(ReflItemList *items, signed int h, signed int k, + signed int l, const char *mero, signed int *hu, + signed int *ku, signed int *lu); + #endif /* SYMMETRY_H */ |