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

Simplify compare_hkl and check_hkl, remove second to last use of "list types"

parent 1e6a810a
Hide whitespace changes
Inline Side-by-side
src/check_hkl.c
View file @ 976170a5
......@@ -54,7 +54,6 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
double num[NBINS];
int cts[NBINS];
int possible[NBINS];
unsigned int *counted;
unsigned int measurements[NBINS];
unsigned int measured[NBINS];
double total_vol, vol_per_shell;
......@@ -73,6 +72,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
int nout = 0;
Reflection *refl;
RefListIterator *iter;
RefList *counted;
int hmax, kmax, lmax;
double asx, asy, asz;
double bsx, bsy, bsz;
......@@ -100,24 +100,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
mean[i] = 0;
}
/* Iterate over all common reflections and calculate min and max
* resolution */
rmin = +INFINITY; rmax = 0.0;
for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
signed int h, k, l;
double d;
get_indices(refl, &h, &k, &l);
d = resolution(cell, h, k, l) * 2.0;
if ( d > rmax ) rmax = d;
if ( d < rmin ) rmin = d;
}
resolution_limits(list, cell, &rmin, &rmax);
STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9);
/* Widen the range just a little bit */
......@@ -155,7 +138,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
rmins[NBINS-1]/1e9, rmaxs[NBINS-1]/1e9);
/* Count the number of reflections possible in each shell */
counted = new_list_count();
counted = reflist_new();
cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
......@@ -170,7 +153,7 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
signed int hs, ks, ls;
int bin;
d = resolution(cell, h, k, l) * 2.0;
d = 2.0 * resolution(cell, h, k, l);
bin = -1;
for ( i=0; i<NBINS; i++ ) {
......@@ -182,21 +165,21 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
if ( bin == -1 ) continue;
get_asymm(sym, h, k, l, &hs, &ks, &ls);
if ( lookup_count(counted, hs, ks, ls) ) continue;
set_count(counted, hs, ks, ls, 1);
if ( find_refl(counted, hs, ks, ls) != NULL ) continue;
add_refl(counted, hs, ks, ls);
possible[bin]++;
}
}
}
free(counted);
reflist_free(counted);
/* Calculate means */
for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
refl = next_refl(refl, iter) )
{
signed int h, k, l;
double d;
int bin;
......@@ -213,14 +196,10 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
break;
}
}
if ( bin == -1 ) {
nout++;
continue;
}
if ( bin == -1 ) continue;
measured[bin]++;
mean[bin] += get_intensity(refl);
}
for ( i=0; i<NBINS; i++ ) {
......@@ -230,8 +209,8 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
/* Characterise the data set */
for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) ) {
refl = next_refl(refl, iter) )
{
signed int h, k, l;
double d;
int bin;
......@@ -266,7 +245,6 @@ static void plot_shells(RefList *list, UnitCell *cell, const SymOpList *sym,
nmeastot += get_redundancy(refl);
var[bin] += pow(val-mean[bin], 2.0);
}
STATUS("overall <snr> = %f\n", snr_total/(double)nmeas);
STATUS("%i measurements in total.\n", nmeastot);
......
src/compare_hkl.c
View file @ 976170a5
......@@ -50,33 +50,24 @@ static void show_help(const char *s)
}
static void plot_shells(RefList *list1, double *arr2, double scale,
UnitCell *cell, const SymOpList *sym,
double rmin_fix, double rmax_fix)
static void plot_shells(RefList *list1, RefList *list2, double scale,
UnitCell *cell, double rmin_fix, double rmax_fix)
{
double num[NBINS];
int cts[NBINS];
int possible[NBINS];
unsigned int *counted;
unsigned int measurements[NBINS];
unsigned int measured[NBINS];
double total_vol, vol_per_shell;
double rmins[NBINS];
double rmaxs[NBINS];
double snr[NBINS];
double den;
double rmin, rmax;
signed int h, k, l;
int i;
int ctot;
FILE *fh;
int nout = 0;
Reflection *refl1;
RefListIterator *iter;
double asx, asy, asz;
double bsx, bsy, bsz;
double csx, csy, csz;
int hmax, kmax, lmax;
FILE *fh;
double den;
int ctot, nout;
if ( cell == NULL ) {
ERROR("Need the unit cell to plot resolution shells.\n");
......@@ -86,30 +77,13 @@ static void plot_shells(RefList *list1, double *arr2, double scale,
for ( i=0; i<NBINS; i++ ) {
num[i] = 0.0;
cts[i] = 0;
possible[i] = 0;
measured[i] = 0;
measurements[i] = 0;
snr[i] = 0;
}
/* Iterate over all common reflections and calculate min and max
* resolution */
rmin = +INFINITY; rmax = 0.0;
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) ) {
signed int h, k, l;
double d;
get_indices(refl1, &h, &k, &l);
d = resolution(cell, h, k, l) * 2.0;
if ( d > rmax ) rmax = d;
if ( d < rmin ) rmin = d;
}
/* Find resolution limits */
resolution_limits(list1, cell, &rmin, &rmax);
STATUS("1/d goes from %f to %f nm^-1\n", rmin/1e9, rmax/1e9);
/* Widen the range just a little bit */
......@@ -147,62 +121,20 @@ static void plot_shells(RefList *list1, double *arr2, double scale,
STATUS("Shell %i: %f to %f\n", NBINS-1,
rmins[NBINS-1]/1e9, rmaxs[NBINS-1]/1e9);
/* Count the number of reflections possible in each shell */
counted = new_list_count();
cell_get_reciprocal(cell, &asx, &asy, &asz,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
hmax = rmax / modulus(asx, asy, asz);
kmax = rmax / modulus(bsx, bsy, bsz);
lmax = rmax / modulus(csx, csy, csz);
for ( h=-hmax; h<hmax; h++ ) {
for ( k=-kmax; k<kmax; k++ ) {
for ( l=-lmax; l<lmax; l++ ) {
double d;
signed int hs, ks, ls;
int bin;
get_asymm(sym, h, k, l, &hs, &ks, &ls);
if ( lookup_count(counted, hs, ks, ls) ) continue;
set_count(counted, hs, ks, ls, 1);
d = resolution(cell, h, k, l) * 2.0;
bin = -1;
for ( i=0; i<NBINS; i++ ) {
if ( (d>rmins[i]) && (d<=rmaxs[i]) ) {
bin = i;
break;
}
}
if ( bin == -1 ) continue;
possible[bin]++;
}
}
}
free(counted);
den = 0.0;
ctot = 0;
nout = 0;
den = 0.0; ctot = 0; nout = 0;
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) ) {
refl1 = next_refl(refl1, iter) )
{
signed int h, k, l;
double d;
int bin;
double i1, i2;
int j;
Reflection *refl2;
get_indices(refl1, &h, &k, &l);
d = resolution(cell, h, k, l) * 2.0;
d = 2.0 * resolution(cell, h, k, l);
bin = -1;
for ( j=0; j<NBINS; j++ ) {
......@@ -218,14 +150,16 @@ static void plot_shells(RefList *list1, double *arr2, double scale,
continue;
}
refl2 = find_refl(list2, h, k, l);
if ( refl2 == NULL ) continue;
i1 = get_intensity(refl1);
i2 = scale * lookup_intensity(arr2, h, k, l);
i2 = scale * get_intensity(refl2);
num[bin] += fabs(i1 - i2);
den += i1 + i2;
ctot++;
cts[bin]++;
}
if ( nout ) {
......@@ -270,18 +204,15 @@ int main(int argc, char *argv[])
int ncom;
RefList *list1;
RefList *list2;
RefList *list1_transformed;
RefList *list2_transformed;
RefList *list1_raw;
RefList *list2_raw;
RefList *ratio;
int config_shells = 0;
char *pdb = NULL;
int rej1 = 0;
int rej2 = 0;
float rmin_fix = -1.0;
float rmax_fix = -1.0;
Reflection *refl1;
RefListIterator *iter;
double *arr2;
int config_unity = 0;
/* Long options */
......@@ -366,88 +297,59 @@ int main(int argc, char *argv[])
cell = load_cell_from_pdb(pdb);
free(pdb);
list1 = read_reflections(afile);
if ( list1 == NULL ) {
list1_raw = read_reflections(afile);
if ( list1_raw == NULL ) {
ERROR("Couldn't read file '%s'\n", afile);
return 1;
}
list2 = read_reflections(bfile);
if ( list2 == NULL ) {
list2_raw = read_reflections(bfile);
if ( list2_raw == NULL ) {
ERROR("Couldn't read file '%s'\n", bfile);
return 1;
}
/* Check that the intensities have the correct symmetry */
if ( check_list_symmetry(list1, sym) ) {
if ( check_list_symmetry(list1_raw, sym) ) {
ERROR("The first input reflection list does not appear to"
" have symmetry %s\n", symmetry_name(sym));
return 1;
}
if ( check_list_symmetry(list2, sym) ) {
if ( check_list_symmetry(list2_raw, sym) ) {
ERROR("The second input reflection list does not appear to"
" have symmetry %s\n", symmetry_name(sym));
return 1;
}
/* Find common reflections (taking symmetry into account) */
list1_transformed = reflist_new();
list2_transformed = reflist_new();
list1 = asymmetric_indices(list1_raw, sym);
list2 = asymmetric_indices(list2_raw, sym);
/* Find common reflections and calculate ratio */
ratio = reflist_new();
ncom = 0;
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) )
{
signed int h, k, l;
signed int he, ke, le;
double val1, val2;
double sig1, sig2;
int ig = 0;
double d;
Reflection *refl2;
Reflection *tr;
get_indices(refl1, &h, &k, &l);
if ( !find_equiv_in_list(list2, h, k, l, sym, &he, &ke, &le) ) {
continue;
}
refl2 = find_refl(list2, h, k, l);
if ( refl2 == NULL ) continue;
copy_data(add_refl(list1_transformed, h, k, l), refl1);
refl2 = find_refl(list2, he, ke, le);
assert(refl2 != NULL);
ncom++;
val1 = get_intensity(refl1);
val2 = get_intensity(refl2);
sig1 = get_esd_intensity(refl1);
sig2 = get_esd_intensity(refl2);
if ( val1 < 3.0 * sig1 ) {
rej1++;
ig = 1;
}
if ( val2 < 3.0 * sig2 ) {
rej2++;
ig = 1;
}
d = 0.5/resolution(cell, h, k, l);
if ( d > 55.0e-10 ) ig = 1;
//if ( d < 15.0e-10 ) ig = 1;
//if ( ig ) continue;
/* Add the old data from 'refl2' to a new list with the same
* indices as its equivalent in 'list1' */
tr = add_refl(list2_transformed, h, k, l);
copy_data(tr, refl2);
/* Add divided version to 'output' list */
tr = add_refl(ratio, h, k, l);
set_int(tr, val1/val2);
set_redundancy(tr, 1);
}
if ( ratiofile != NULL ) {
......@@ -457,67 +359,55 @@ int main(int argc, char *argv[])
gsl_set_error_handler_off();
STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej1, afile);
STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej2, bfile);
ncom = num_reflections(list2_transformed);
STATUS("%i,%i reflections: %i in common\n",
num_reflections(list1), num_reflections(list2), ncom);
reflist_free(list1);
reflist_free(list2);
/* Trimming the other way round is not necessary,
* because of these two lines */
arr2 = intensities_from_list(list2_transformed);
reflist_free(list2_transformed);
R1 = stat_r1_ignore(list1_transformed, arr2, &scale_r1fi, config_unity);
R1 = stat_r1_ignore(list1, list2, &scale_r1fi, config_unity);
STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
" (ignoring negative intensities)\n",
R1*100.0, scale_r1fi);
R1 = stat_r1_zero(list1_transformed, arr2, &scale_r1, config_unity);
R1 = stat_r1_zero(list1, list2, &scale_r1, config_unity);
STATUS("R1(F) = %5.4f %% (scale=%5.2e)"
" (zeroing negative intensities)\n",
R1*100.0, scale_r1);
R2 = stat_r2(list1_transformed, arr2, &scale_r2, config_unity);
R2 = stat_r2(list1, list2, &scale_r2, config_unity);
STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
R1i = stat_r1_i(list1_transformed, arr2, &scale_r1i, config_unity);
R1i = stat_r1_i(list1, list2, &scale_r1i, config_unity);
STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i);
Rdiff = stat_rdiff_ignore(list1_transformed, arr2, &scale_rdig,
Rdiff = stat_rdiff_ignore(list1, list2, &scale_rdig,
config_unity);
STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
" (ignoring negative intensities)\n",
Rdiff*100.0, scale_rdig);
Rdiff = stat_rdiff_zero(list1_transformed, arr2, &scale, config_unity);
Rdiff = stat_rdiff_zero(list1, list2, &scale, config_unity);
STATUS("Rint(F) = %5.4f %% (scale=%5.2e)"
" (zeroing negative intensities)\n",
Rdiff*100.0, scale);
Rdiff = stat_rdiff_intensity(list1_transformed, arr2, &scale_rintint,
Rdiff = stat_rdiff_intensity(list1, list2, &scale_rintint,
config_unity);
STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n",
Rdiff*100.0, scale_rintint);
pearson = stat_pearson_i(list1_transformed, arr2);
pearson = stat_pearson_i(list1, list2);
STATUS("Pearson r(I) = %5.4f\n", pearson);
pearson = stat_pearson_f_ignore(list1_transformed, arr2);
pearson = stat_pearson_f_ignore(list1, list2);
STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
pearson);
pearson = stat_pearson_f_zero(list1_transformed, arr2);
pearson = stat_pearson_f_zero(list1, list2);
STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
pearson);
if ( config_shells ) {
plot_shells(list1_transformed, arr2, scale_r1i,
cell, sym, rmin_fix, rmax_fix);
plot_shells(list1, list2, scale_r1i,
cell, rmin_fix, rmax_fix);
}
return 0;
......
src/statistics.c
View file @ 976170a5
......@@ -38,7 +38,7 @@
struct r_params {
RefList *list1;
double *arr2;
RefList *list2;
int fom; /* Which FoM to use (see the enum just below) */
};
......@@ -54,9 +54,9 @@ enum {
/* Return the least squares optimal scaling factor when comparing intensities.
* list1,arr2 are the two intensity lists to compare.
* list1,list2 are the two intensity lists to compare.
*/
double stat_scale_intensity(RefList *list1, double *arr2)
double stat_scale_intensity(RefList *list1, RefList *list2)
{
double top = 0.0;
double bot = 0.0;
......@@ -65,15 +65,18 @@ double stat_scale_intensity(RefList *list1, double *arr2)
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) ) {
refl1 = next_refl(refl1, iter) )
{
double i1, i2;
signed int h, k, l;
Reflection *refl2;
get_indices(refl1, &h, &k, &l);
i2 = lookup_intensity(arr2, h, k, l);
refl2 = find_refl(list2, h, k, l);
if ( refl2 == NULL ) continue; /* No common reflection */
i1 = get_intensity(refl1);
i2 = get_intensity(refl2);
top += i1 * i2;
bot += i2 * i2;
......@@ -87,10 +90,10 @@ double stat_scale_intensity(RefList *list1, double *arr2)
/* Return the least squares optimal scaling factor when comparing the square
* roots of the intensities (i.e. one approximation to the structure factor
* moduli).
* list1,arr2 are the two intensity lists to compare (they contain intensities,
* list1,list2 are the two intensity lists to compare (they contain intensities,
* not square rooted intensities).
*/
static double stat_scale_sqrti(RefList *list1, double *arr2)
static double stat_scale_sqrti(RefList *list1, RefList *list2)
{
double top = 0.0;
double bot = 0.0;
......@@ -99,16 +102,19 @@ static double stat_scale_sqrti(RefList *list1, double *arr2)
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) ) {
refl1 = next_refl(refl1, iter) )
{
double i1, i2;
double f1, f2;
signed int h, k, l;
Reflection *refl2;
get_indices(refl1, &h, &k, &l);
i2 = lookup_intensity(arr2, h, k, l);
refl2 = find_refl(list2, h, k, l);
if ( refl2 == NULL ) continue; /* No common reflection */
i1 = get_intensity(refl1);
i2 = get_intensity(refl2);
if ( i1 < 0.0 ) continue;
f1 = sqrt(i1);
......@@ -125,7 +131,7 @@ static double stat_scale_sqrti(RefList *list1, double *arr2)
}
static double internal_r1_ignorenegs(RefList *list1, double *arr2,
static double internal_r1_ignorenegs(RefList *list1, RefList *list2,
double scale)
{
double top = 0.0;
......@@ -135,16 +141,19 @@ static double internal_r1_ignorenegs(RefList *list1, double *arr2,
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) ) {
refl1 = next_refl(refl1, iter) )
{
double i1, i2;
double f1, f2;
signed int h, k, l;
Reflection *refl2;
get_indices(refl1, &h, &k, &l);
i2 = lookup_intensity(arr2, h, k, l);
refl2 = find_refl(list2, h, k, l);
if ( refl2 == NULL ) continue; /* No common reflection */
i1 = get_intensity(refl1);
i2 = get_intensity(refl2);
if ( i1 < 0.0 ) continue;
f1 = sqrt(i1);
......@@ -162,7 +171,7 @@ static double internal_r1_ignorenegs(RefList *list1, double *arr2,
}
static double internal_r1_negstozero(RefList *list1, double *arr2,
static double internal_r1_negstozero(RefList *list1, RefList *list2,
double scale)
{
double top = 0.0;
......@@ -172,16 +181,19 @@ static double internal_r1_negstozero(RefList *list1, double *arr2,
for ( refl1 = first_refl(list1, &iter);
refl1 != NULL;
refl1 = next_refl(refl1, iter) ) {
refl1 = next_refl(refl1, iter) )
{
double i1, i2;
double f1, f2;
signed int h, k, l;
Reflection *refl2;
get_indices(refl1, &h, &k, &l);
i2 = lookup_intensity(arr2, h, k, l);
refl2 = find_refl(list2, h, k, l);
if ( refl2 == NULL ) continue; /* No common reflection */
i1 = get_intensity(refl1);
i2 = get_intensity(refl2);
f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
......@@ -197,7 +209,7 @@ static double internal_r1_negstozero(RefList *list1, double *arr2,
}
static double internal_r2(RefList *list1, double *arr2, double scale)