From 976170a5a1838077bfb230dfa634ed8310543815 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 20 Sep 2011 16:56:43 +0200 Subject: Simplify compare_hkl and check_hkl, remove second to last use of "list types" --- src/check_hkl.c | 46 ++++--------- src/compare_hkl.c | 198 ++++++++++++----------------------------------------- src/statistics.c | 201 +++++++++++++++++++++++++++++++----------------------- src/statistics.h | 23 ++++--- 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 = %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 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; hrmins[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 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 */ -- cgit v1.2.3