Commit 003a1a19 authored by Thomas White's avatar Thomas White Committed by Thomas White
Browse files

compare_hkl: Take symmetry into account

parent 84a7190c
......@@ -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 \
......
......@@ -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 \
......
......@@ -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 ) {
......
......@@ -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) ) {
......
......@@ -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;
}
......@@ -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 */
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment