From fbd0ad07cf2bcdc7d48e6b5bc0f4b37c9f7e9817 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 24 Jun 2011 11:40:20 +0200 Subject: Add "-u" (unity scaling) option to compare_hkl --- src/compare_hkl.c | 23 +++++++---- src/statistics.c | 117 +++++++++++++++++++++++++++++------------------------- src/statistics.h | 17 ++++---- 3 files changed, 88 insertions(+), 69 deletions(-) diff --git a/src/compare_hkl.c b/src/compare_hkl.c index 4ab2dc70..02a85d2f 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -277,6 +277,7 @@ int main(int argc, char *argv[]) Reflection *refl1; RefListIterator *iter; double *arr2; + int config_unity = 0; /* Long options */ const struct option longopts[] = { @@ -291,7 +292,9 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "ho:y:p:", longopts, NULL)) != -1) { + while ((c = getopt_long(argc, argv, "ho:y:p:u", + longopts, NULL)) != -1) + { switch (c) { case 'h' : @@ -310,6 +313,10 @@ int main(int argc, char *argv[]) pdb = strdup(optarg); break; + case 'u' : + config_unity = 1; + break; + case 0 : break; @@ -473,33 +480,33 @@ int main(int argc, char *argv[]) arr2 = intensities_from_list(list2_transformed); reflist_free(list2_transformed); - R1 = stat_r1_ignore(list1, arr2, &scale_r1fi); + R1 = stat_r1_ignore(list1, arr2, &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, arr2, &scale_r1); + R1 = stat_r1_zero(list1, arr2, &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, arr2, &scale_r2); + R2 = stat_r2(list1, arr2, &scale_r2, config_unity); STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2); - R1i = stat_r1_i(list1, arr2, &scale_r1i); + R1i = stat_r1_i(list1, arr2, &scale_r1i, config_unity); STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i); - Rdiff = stat_rdiff_ignore(list1, arr2, &scale_rdig); + Rdiff = stat_rdiff_ignore(list1, arr2, &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, arr2, &scale); + Rdiff = stat_rdiff_zero(list1, arr2, &scale, config_unity); STATUS("Rint(F) = %5.4f %% (scale=%5.2e)" " (zeroing negative intensities)\n", Rdiff*100.0, scale); - Rdiff = stat_rdiff_intensity(list1, arr2, &scale_rintint); + Rdiff = stat_rdiff_intensity(list1, arr2, &scale_rintint, config_unity); STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n", Rdiff*100.0, scale_rintint); diff --git a/src/statistics.c b/src/statistics.c index 2312c4a7..791952b8 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -376,7 +376,8 @@ 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, double *arr2, double *scalep, int fom, + int u) { gsl_function F; gsl_min_fminimizer *s; @@ -389,57 +390,65 @@ static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom) rp.arr2 = arr2; rp.fom = fom; - F.function = &calc_r; - F.params = &rp; + if ( u ) { - s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent); + scale = 1.0; - /* Initial guess */ - switch ( fom ) { - case R_1_ZERO : - case R_1_IGNORE : - case R_DIFF_ZERO : - case R_DIFF_IGNORE : - scale = stat_scale_sqrti(list1, arr2); - break; - case R_2 : - case R_1_I : - case R_DIFF_INTENSITY : - scale = stat_scale_intensity(list1, arr2); - break; - } - //STATUS("Initial scale factor estimate: %5.2e\n", scale); - - /* Probably within an order of magnitude either side */ - gsl_min_fminimizer_set(s, &F, scale, scale/10.0, scale*10.0); + } else { - do { + F.function = &calc_r; + F.params = &rp; - double lo, up; + s = gsl_min_fminimizer_alloc(gsl_min_fminimizer_brent); - /* Iterate */ - if ( gsl_min_fminimizer_iterate(s) ) { - ERROR("Failed to find scale factor.\n"); - return NAN; + /* Initial guess */ + switch ( fom ) { + case R_1_ZERO : + case R_1_IGNORE : + case R_DIFF_ZERO : + case R_DIFF_IGNORE : + scale = stat_scale_sqrti(list1, arr2); + break; + case R_2 : + case R_1_I : + case R_DIFF_INTENSITY : + scale = stat_scale_intensity(list1, arr2); + break; } + //STATUS("Initial scale factor estimate: %5.2e\n", scale); - /* Get the current estimate */ - scale = gsl_min_fminimizer_x_minimum(s); - lo = gsl_min_fminimizer_x_lower(s); - up = gsl_min_fminimizer_x_upper(s); + /* Probably within an order of magnitude either side */ + gsl_min_fminimizer_set(s, &F, scale, scale/10.0, scale*10.0); - /* Check for convergence */ - status = gsl_min_test_interval(lo, up, 0.001, 0.0); + do { - iter++; + double lo, up; - } while ( status == GSL_CONTINUE ); + /* Iterate */ + if ( gsl_min_fminimizer_iterate(s) ) { + ERROR("Failed to find scale factor.\n"); + return NAN; + } - if (status != GSL_SUCCESS) { - ERROR("Scale factor minimisation failed.\n"); - } + /* Get the current estimate */ + scale = gsl_min_fminimizer_x_minimum(s); + lo = gsl_min_fminimizer_x_lower(s); + up = gsl_min_fminimizer_x_upper(s); + + /* Check for convergence */ + status = gsl_min_test_interval(lo, up, 0.001, 0.0); + + iter++; + + } while ( status == GSL_CONTINUE ); - gsl_min_fminimizer_free(s); + if ( status != GSL_SUCCESS ) { + ERROR("Scale factor minimisation failed.\n"); + } + + gsl_min_fminimizer_free(s); + + } //STATUS("Final scale factor: %5.2e\n", scale); *scalep = scale; @@ -447,45 +456,45 @@ static double r_minimised(RefList *list1, double *arr2, double *scalep, int fom) } -double stat_r1_ignore(RefList *list1, double *arr2, double *scalep) +double stat_r1_ignore(RefList *list1, double *arr2, double *scalep, int u) { - return r_minimised(list1, arr2, scalep, R_1_IGNORE); + return r_minimised(list1, arr2, scalep, R_1_IGNORE, u); } -double stat_r1_zero(RefList *list1, double *arr2, double *scalep) +double stat_r1_zero(RefList *list1, double *arr2, double *scalep, int u) { - return r_minimised(list1, arr2, scalep, R_1_ZERO); + return r_minimised(list1, arr2, scalep, R_1_ZERO, u); } -double stat_r2(RefList *list1, double *arr2, double *scalep) +double stat_r2(RefList *list1, double *arr2, double *scalep, int u) { - return r_minimised(list1, arr2, scalep, R_2); + return r_minimised(list1, arr2, scalep, R_2, u); } -double stat_r1_i(RefList *list1, double *arr2, double *scalep) +double stat_r1_i(RefList *list1, double *arr2, double *scalep, int u) { - return r_minimised(list1, arr2, scalep, R_1_I); + return r_minimised(list1, arr2, scalep, R_1_I, u); } -double stat_rdiff_zero(RefList *list1, double *arr2, double *scalep) +double stat_rdiff_zero(RefList *list1, double *arr2, double *scalep, int u) { - return r_minimised(list1, arr2, scalep, R_DIFF_ZERO); + return r_minimised(list1, arr2, scalep, R_DIFF_ZERO, u); } -double stat_rdiff_ignore(RefList *list1, double *arr2, double *scalep) +double stat_rdiff_ignore(RefList *list1, double *arr2, double *scalep, int u) { - return r_minimised(list1, arr2, scalep, R_DIFF_IGNORE); + return r_minimised(list1, arr2, scalep, R_DIFF_IGNORE, u); } -double stat_rdiff_intensity(RefList *list1, double *arr2, double *scalep) +double stat_rdiff_intensity(RefList *list1, double *arr2, double *scalep, int u) { - return r_minimised(list1, arr2, scalep, R_DIFF_INTENSITY); + return r_minimised(list1, arr2, scalep, R_DIFF_INTENSITY, u); } diff --git a/src/statistics.h b/src/statistics.h index e1d7244e..ced295de 100644 --- a/src/statistics.h +++ b/src/statistics.h @@ -21,17 +21,20 @@ extern double stat_scale_intensity(RefList *list1, double *arr2); -extern double stat_r1_zero(RefList *list1, double *arr2, double *scalep); -extern double stat_r1_ignore(RefList *list1, double *arr2, double *scalep); +extern double stat_r1_zero(RefList *list1, double *arr2, double *scalep, int u); +extern double stat_r1_ignore(RefList *list1, double *arr2, + double *scalep, int u); -extern double stat_r2(RefList *list1, double *arr2, double *scalep); +extern double stat_r2(RefList *list1, double *arr2, double *scalep, int u); -extern double stat_r1_i(RefList *list1, double *arr2, double *scalep); +extern double stat_r1_i(RefList *list1, double *arr2, double *scalep, int u); -extern double stat_rdiff_zero(RefList *list1, double *arr2, double *scalep); -extern double stat_rdiff_ignore(RefList *list1, double *arr2, double *scalep); +extern double stat_rdiff_zero(RefList *list1, double *arr2, + double *scalep, int u); +extern double stat_rdiff_ignore(RefList *list1, double *arr2, + double *scalep, int u); extern double stat_rdiff_intensity(RefList *list1, double *arr2, - double *scalep); + double *scalep, int u); extern double stat_pearson_i(RefList *list1, double *arr2); extern double stat_pearson_f_zero(RefList *list1, double *arr2); -- cgit v1.2.3