aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-09-20 16:56:43 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:37 +0100
commit976170a5a1838077bfb230dfa634ed8310543815 (patch)
treea55ee0ec741b7604a1abf5f5f0274f7ffc8ff176 /src
parent1e6a810ad46056154dd8984773d15828838658b5 (diff)
Simplify compare_hkl and check_hkl, remove second to last use of "list types"
Diffstat (limited to 'src')
-rw-r--r--src/check_hkl.c46
-rw-r--r--src/compare_hkl.c198
-rw-r--r--src/statistics.c201
-rw-r--r--src/statistics.h23
4 files changed, 185 insertions, 283 deletions
diff --git a/src/check_hkl.c b/src/check_hkl.c
index 6eb1d26e..c9d50c20 100644
--- a/src/check_hkl.c
+++ b/src/check_hkl.c
@@ -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);
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index c61df8c4..7df100c0 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -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;
diff --git a/src/statistics.c b/src/statistics.c
index d9652ece..7990672d 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -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)
+static double internal_r2(RefList *list1, RefList *list2, double scale)
{
double top = 0.0;
double bot = 0.0;
@@ -206,15 +218,19 @@ static double internal_r2(RefList *list1, double *arr2, double scale)
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);
+
i2 *= scale;
top += pow(i1 - i2, 2.0);
@@ -226,7 +242,7 @@ static double internal_r2(RefList *list1, double *arr2, double scale)
}
-static double internal_r_i(RefList *list1, double *arr2, double scale)
+static double internal_r_i(RefList *list1, RefList *list2, double scale)
{
double top = 0.0;
double bot = 0.0;
@@ -235,15 +251,18 @@ static double internal_r_i(RefList *list1, double *arr2, double scale)
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);
i2 *= scale;
top += fabs(i1-i2);
@@ -255,7 +274,7 @@ static double internal_r_i(RefList *list1, double *arr2, double scale)
}
-static double internal_rdiff_intensity(RefList *list1, double *arr2,
+static double internal_rdiff_intensity(RefList *list1, RefList *list2,
double scale)
{
double top = 0.0;
@@ -265,16 +284,18 @@ static double internal_rdiff_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);
i2 *= scale;
top += fabs(i1 - i2);
@@ -286,7 +307,7 @@ static double internal_rdiff_intensity(RefList *list1, double *arr2,
}
-static double internal_rdiff_negstozero(RefList *list1, double *arr2,
+static double internal_rdiff_negstozero(RefList *list1, RefList *list2,
double scale)
{
double top = 0.0;
@@ -296,16 +317,19 @@ static double internal_rdiff_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;
@@ -321,7 +345,7 @@ static double internal_rdiff_negstozero(RefList *list1, double *arr2,
}
-static double internal_rdiff_ignorenegs(RefList *list1, double *arr2,
+static double internal_rdiff_ignorenegs(RefList *list1, RefList *list2,
double scale)
{
double top = 0.0;
@@ -331,16 +355,19 @@ static double internal_rdiff_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);
@@ -362,23 +389,23 @@ static double calc_r(double scale, void *params)
{
struct r_params *rp = params;
- switch ( rp->fom) {
+ switch ( rp->fom ) {
case R_1_ZERO :
- return internal_r1_negstozero(rp->list1, rp->arr2, scale);
+ return internal_r1_negstozero(rp->list1, rp->list2, scale);
case R_1_IGNORE :
- return internal_r1_ignorenegs(rp->list1, rp->arr2, scale);
+ return internal_r1_ignorenegs(rp->list1, rp->list2, scale);
case R_2 :
- return internal_r2(rp->list1, rp->arr2, scale);
+ return internal_r2(rp->list1, rp->list2, scale);
case R_1_I :
- return internal_r_i(rp->list1, rp->arr2, scale);
+ return internal_r_i(rp->list1, rp->list2, scale);
case R_DIFF_ZERO :
- return internal_rdiff_negstozero(rp->list1, rp->arr2,scale);
+ return internal_rdiff_negstozero(rp->list1, rp->list2,scale);
case R_DIFF_IGNORE :
- return internal_rdiff_ignorenegs(rp->list1, rp->arr2, scale);
+ return internal_rdiff_ignorenegs(rp->list1, rp->list2, scale);
case R_DIFF_INTENSITY :
- return internal_rdiff_intensity(rp->list1, rp->arr2, scale);
+ return internal_rdiff_intensity(rp->list1, rp->list2, scale);
}
ERROR("No such FoM!\n");
@@ -386,7 +413,7 @@ static double calc_r(double scale, void *params)
}
-static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom,
+static double r_minimised(RefList *list1, RefList *list2, double *scalep, int fom,
int u)
{
gsl_function F;
@@ -397,7 +424,7 @@ static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom,
int iter = 0;
rp.list1 = list1;
- rp.arr2 = arr2;
+ rp.list2 = list2;
rp.fom = fom;
if ( u ) {
@@ -417,12 +444,12 @@ static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom,
case R_1_IGNORE :
case R_DIFF_ZERO :
case R_DIFF_IGNORE :
- scale = stat_scale_sqrti(list1, arr2);
+ scale = stat_scale_sqrti(list1, list2);
break;
case R_2 :
case R_1_I :
case R_DIFF_INTENSITY :
- scale = stat_scale_intensity(list1, arr2);
+ scale = stat_scale_intensity(list1, list2);
break;
}
//STATUS("Initial scale factor estimate: %5.2e\n", scale);
@@ -466,49 +493,49 @@ static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom,
}
-double stat_r1_ignore(RefList *list1, double *arr2, double *scalep, int u)
+double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep, int u)
{
- return r_minimised(list1, arr2, scalep, R_1_IGNORE, u);
+ return r_minimised(list1, list2, scalep, R_1_IGNORE, u);
}
-double stat_r1_zero(RefList *list1, double *arr2, double *scalep, int u)
+double stat_r1_zero(RefList *list1, RefList *list2, double *scalep, int u)
{
- return r_minimised(list1, arr2, scalep, R_1_ZERO, u);
+ return r_minimised(list1, list2, scalep, R_1_ZERO, u);
}
-double stat_r2(RefList *list1, double *arr2, double *scalep, int u)
+double stat_r2(RefList *list1, RefList *list2, double *scalep, int u)
{
- return r_minimised(list1, arr2, scalep, R_2, u);
+ return r_minimised(list1, list2, scalep, R_2, u);
}
-double stat_r1_i(RefList *list1, double *arr2, double *scalep, int u)
+double stat_r1_i(RefList *list1, RefList *list2, double *scalep, int u)
{
- return r_minimised(list1, arr2, scalep, R_1_I, u);
+ return r_minimised(list1, list2, scalep, R_1_I, u);
}
-double stat_rdiff_zero(RefList *list1, double *arr2, double *scalep, int u)
+double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep, int u)
{
- return r_minimised(list1, arr2, scalep, R_DIFF_ZERO, u);
+ return r_minimised(list1, list2, scalep, R_DIFF_ZERO, u);
}
-double stat_rdiff_ignore(RefList *list1, double *arr2, double *scalep, int u)
+double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep, int u)
{
- return r_minimised(list1, arr2, scalep, R_DIFF_IGNORE, u);
+ return r_minimised(list1, list2, scalep, R_DIFF_IGNORE, u);
}
-double stat_rdiff_intensity(RefList *list1, double *arr2, double *scalep, int u)
+double stat_rdiff_intensity(RefList *list1, RefList *list2, double *scalep, int u)
{
- return r_minimised(list1, arr2, scalep, R_DIFF_INTENSITY, u);
+ return r_minimised(list1, list2, scalep, R_DIFF_INTENSITY, u);
}
-double stat_pearson_i(RefList *list1, double *arr2)
+double stat_pearson_i(RefList *list1, RefList *list2)
{
double *vec1, *vec2;
int ni = num_reflections(list1);
@@ -522,15 +549,18 @@ double stat_pearson_i(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);
vec1[nacc] = i1;
vec2[nacc] = i2;
@@ -545,7 +575,7 @@ double stat_pearson_i(RefList *list1, double *arr2)
}
-double stat_pearson_f_ignore(RefList *list1, double *arr2)
+double stat_pearson_f_ignore(RefList *list1, RefList *list2)
{
double *vec1, *vec2;
int ni = num_reflections(list1);
@@ -559,26 +589,26 @@ double stat_pearson_f_ignore(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;
- int ig = 0;
+ 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 ) ig = 1;
- f1 = sqrt(i1);
+ if ( i1 < 0.0 ) continue;
+ if ( i2 < 0.0 ) continue;
- if ( i2 < 0.0 ) ig = 1;
+ f1 = sqrt(i1);
f2 = sqrt(i2);
- if ( ig ) continue;
-
vec1[nacc] = f1;
vec2[nacc] = f2;
nacc++;
@@ -593,7 +623,7 @@ double stat_pearson_f_ignore(RefList *list1, double *arr2)
}
-double stat_pearson_f_zero(RefList *list1, double *arr2)
+double stat_pearson_f_zero(RefList *list1, RefList *list2)
{
double *vec1, *vec2;
int ni = num_reflections(list1);
@@ -607,16 +637,19 @@ double stat_pearson_f_zero(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;
f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
diff --git a/src/statistics.h b/src/statistics.h
index ced295de..8fa78ea6 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -19,26 +19,27 @@
#include "reflist.h"
-extern double stat_scale_intensity(RefList *list1, double *arr2);
+extern double stat_scale_intensity(RefList *list1, RefList *list2);
-extern double stat_r1_zero(RefList *list1, double *arr2, double *scalep, int u);
-extern double stat_r1_ignore(RefList *list1, double *arr2,
+extern double stat_r1_zero(RefList *list1, RefList *list2,
+ double *scalep, int u);
+extern double stat_r1_ignore(RefList *list1, RefList *list2,
double *scalep, int u);
-extern double stat_r2(RefList *list1, double *arr2, double *scalep, int u);
+extern double stat_r2(RefList *list1, RefList *list2, double *scalep, int u);
-extern double stat_r1_i(RefList *list1, double *arr2, double *scalep, int u);
+extern double stat_r1_i(RefList *list1, RefList *list2, double *scalep, int u);
-extern double stat_rdiff_zero(RefList *list1, double *arr2,
+extern double stat_rdiff_zero(RefList *list1, RefList *list2,
double *scalep, int u);
-extern double stat_rdiff_ignore(RefList *list1, double *arr2,
+extern double stat_rdiff_ignore(RefList *list1, RefList *list2,
double *scalep, int u);
-extern double stat_rdiff_intensity(RefList *list1, double *arr2,
+extern double stat_rdiff_intensity(RefList *list1, RefList *list2,
double *scalep, int u);
-extern double stat_pearson_i(RefList *list1, double *arr2);
-extern double stat_pearson_f_zero(RefList *list1, double *arr2);
-extern double stat_pearson_f_ignore(RefList *list1, double *arr2);
+extern double stat_pearson_i(RefList *list1, RefList *list2);
+extern double stat_pearson_f_zero(RefList *list1, RefList *list2);
+extern double stat_pearson_f_ignore(RefList *list1, RefList *list2);
#endif /* STATISTICS_H */