aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/compare_hkl.c23
-rw-r--r--src/statistics.c117
-rw-r--r--src/statistics.h17
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);
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-05 17:40:31 +0000