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

Add "-u" (unity scaling) option to compare_hkl

parent dca4929e
Hide whitespace changes
Inline Side-by-side
src/compare_hkl.c
View file @ fbd0ad07
......@@ -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);
......
src/statistics.c
View file @ fbd0ad07
......@@ -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);
}
......
src/statistics.h
View file @ fbd0ad07
......@@ -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);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

imprint | privacy | Conditions of Use