aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/compare_hkl.c39
-rw-r--r--src/process_hkl.c2
-rw-r--r--src/statistics.c257
-rw-r--r--src/statistics.h21
4 files changed, 283 insertions, 36 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index ead7b062..272b1672 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -137,7 +137,7 @@ int main(int argc, char *argv[])
char *afile = NULL;
char *bfile = NULL;
char *sym = NULL;
- double scale, scale_r2, R1, R2, Rdiff, pearson;
+ double scale, scale_r2, scale_rint, R1, R2, Rint, Rdiff, pearson;
int i, ncom;
ReflItemList *i1, *i2, *icommon;
int config_luzzati = 0;
@@ -245,14 +245,39 @@ int main(int argc, char *argv[])
STATUS("%i,%i reflections: %i in common\n",
num_items(i1), num_items(i2), ncom);
- R1 = stat_r1(ref1, ref2_transformed, icommon, &scale);
- STATUS("R1 = %5.4f %% (scale=%5.2e)\n", R1*100.0, scale);
+
+ R1 = stat_r1_ignore(ref1, ref2_transformed, icommon, &scale);
+ STATUS("R1 = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
+ R1*100.0, scale);
+
+ R1 = stat_r1_zero(ref1, ref2_transformed, icommon, &scale);
+ STATUS("R1 = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
+ R1*100.0, scale);
+
R2 = stat_r2(ref1, ref2_transformed, icommon, &scale_r2);
STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2);
- Rdiff = stat_rdiff(ref1, ref2_transformed, icommon, &scale);
- STATUS("Rdiff = %5.4f %% (scale=%5.2e)\n", Rdiff*100.0, scale);
- pearson = stat_pearson(ref1, ref2_transformed, icommon);
- STATUS("Pearson r = %5.4f\n", pearson);
+
+ Rint = stat_rint(ref1, ref2_transformed, icommon, &scale_rint);
+ STATUS("Rint = %5.4f %% (scale=%5.2e)\n", Rint*100.0, scale_rint);
+
+ Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale);
+ STATUS("Rdiff = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n",
+ Rdiff*100.0, scale);
+
+ Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale);
+ STATUS("Rdiff = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n",
+ Rdiff*100.0, scale);
+
+ pearson = stat_pearson_i(ref1, ref2_transformed, icommon);
+ STATUS("Pearson r(I) = %5.4f\n", pearson);
+
+ pearson = stat_pearson_f_ignore(ref1, ref2_transformed, icommon);
+ STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n",
+ pearson);
+
+ pearson = stat_pearson_f_zero(ref1, ref2_transformed, icommon);
+ STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n",
+ pearson);
if ( config_luzzati ) {
plot_luzzati(ref1, ref2_transformed, icommon, scale_r2, cell);
diff --git a/src/process_hkl.c b/src/process_hkl.c
index beeebb39..644fdfea 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -187,7 +187,7 @@ static int resolve_twin(const double *model, ReflItemList *observed,
}
intersection = intersection_items(observed, items);
- fom = stat_pearson(trial_ints, model, intersection);
+ fom = stat_pearson_i(trial_ints, model, intersection);
delete_items(intersection);
free(trial_ints);
diff --git a/src/statistics.c b/src/statistics.c
index 81b84598..a6936798 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -32,9 +32,12 @@ struct r_params {
};
enum {
- R_1,
+ R_1_ZERO,
+ R_1_IGNORE,
R_2,
- R_DIFF,
+ R_INT,
+ R_DIFF_ZERO,
+ R_DIFF_IGNORE,
};
@@ -94,9 +97,11 @@ static double stat_scale_sqrti(const double *ref1, const double *ref2,
h = it->h; k = it->k; l = it->l;
i1 = lookup_intensity(ref1, h, k, l);
- f1 = i1 > 0.0 ? sqrt(i1) : 0.0;
+ if ( i1 < 0.0 ) continue;
+ f1 = sqrt(i1);
i2 = lookup_intensity(ref2, h, k, l);
- f2 = i2 > 0.0 ? sqrt(i2) : 0.0;
+ if ( i2 < 0.0 ) continue;
+ f2 = sqrt(i2);
top += f1 * f2;
bot += f2 * f2;
@@ -107,8 +112,41 @@ static double stat_scale_sqrti(const double *ref1, const double *ref2,
}
-static double internal_r1(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_r1_ignorenegs(const double *ref1, const double *ref2,
+ ReflItemList *items, double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ int i;
+
+ for ( i=0; i<num_items(items); i++ ) {
+
+ double i1, f1, i2, f2;
+ struct refl_item *it;
+ signed int h, k, l;
+
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
+
+ i1 = lookup_intensity(ref1, h, k, l);
+ if ( i1 < 0.0 ) continue;
+ f1 = sqrt(i1);
+ i2 = lookup_intensity(ref2, h, k, l);
+ if ( i2 < 0.0 ) continue;
+ f2 = sqrt(i2);
+ f2 *= scale;
+
+ top += fabs(f1 - f2);
+ bot += f1;
+
+ }
+
+ return top/bot;
+}
+
+
+static double internal_r1_negstozero(const double *ref1, const double *ref2,
+ ReflItemList *items, double scale)
{
double top = 0.0;
double bot = 0.0;
@@ -166,8 +204,36 @@ static double internal_r2(const double *ref1, const double *ref2,
}
-static double internal_rdiff(const double *ref1, const double *ref2,
- ReflItemList *items, double scale)
+static double internal_rint(const double *ref1, const double *ref2,
+ ReflItemList *items, double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ int i;
+
+ for ( i=0; i<num_items(items); i++ ) {
+
+ double i1, i2;
+ struct refl_item *it;
+ signed int h, k, l;
+
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
+
+ i1 = lookup_intensity(ref1, h, k, l);
+ i2 = scale * lookup_intensity(ref2, h, k, l);
+
+ top += fabs(i1-i2);
+ bot += fabs(i1);
+
+ }
+
+ return top/bot;
+}
+
+
+static double internal_rdiff_negstozero(const double *ref1, const double *ref2,
+ ReflItemList *items, double scale)
{
double top = 0.0;
double bot = 0.0;
@@ -197,18 +263,62 @@ static double internal_rdiff(const double *ref1, const double *ref2,
}
+static double internal_rdiff_ignorenegs(const double *ref1, const double *ref2,
+ ReflItemList *items, double scale)
+{
+ double top = 0.0;
+ double bot = 0.0;
+ int i;
+
+ for ( i=0; i<num_items(items); i++ ) {
+
+ double i1, i2, f1, f2;
+ struct refl_item *it;
+ signed int h, k, l;
+
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
+
+ i1 = lookup_intensity(ref1, h, k, l);
+ if ( i1 < 0.0 ) continue;
+ f1 = sqrt(i1);
+ i2 = lookup_intensity(ref2, h, k, l);
+ if ( i2 < 0.0 ) continue;
+ f2 = sqrt(i2);
+ f2 *= scale;
+
+ top += fabs(f1 - f2);
+ bot += f1 + f2;
+
+ }
+
+ return 2.0*top/bot;
+}
+
static double calc_r(double scale, void *params)
{
struct r_params *rp = params;
switch ( rp->fom) {
- case R_1 :
- return internal_r1(rp->ref1, rp->ref2, rp->items, scale);
+ case R_1_ZERO :
+ return internal_r1_negstozero(rp->ref1, rp->ref2,
+ rp->items, scale);
+ case R_1_IGNORE :
+ return internal_r1_ignorenegs(rp->ref1, rp->ref2,
+ rp->items, scale);
case R_2 :
return internal_r2(rp->ref1, rp->ref2, rp->items, scale);
- case R_DIFF :
- return internal_rdiff(rp->ref1, rp->ref2, rp->items, scale);
+
+ case R_INT :
+ return internal_rint(rp->ref1, rp->ref2, rp->items, scale);
+
+ case R_DIFF_ZERO :
+ return internal_rdiff_negstozero(rp->ref1, rp->ref2,
+ rp->items, scale);
+ case R_DIFF_IGNORE :
+ return internal_rdiff_ignorenegs(rp->ref1, rp->ref2,
+ rp->items, scale);
}
ERROR("No such FoM!\n");
@@ -238,15 +348,16 @@ static double r_minimised(const double *ref1, const double *ref2,
/* Initial guess */
switch ( fom ) {
- case R_1 :
+ case R_1_ZERO :
+ case R_1_IGNORE :
+ case R_DIFF_ZERO :
+ case R_DIFF_IGNORE :
scale = stat_scale_sqrti(ref1, ref2, items);
break;
case R_2 :
+ case R_INT :
scale = stat_scale_intensity(ref1, ref2, items);
break;
- case R_DIFF :
- scale = stat_scale_sqrti(ref1, ref2, items);
- break;
}
//STATUS("Initial scale factor estimate: %5.2e\n", scale);
@@ -284,10 +395,17 @@ static double r_minimised(const double *ref1, const double *ref2,
}
-double stat_r1(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_r1_ignore(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep)
{
- return r_minimised(ref1, ref2, items, scalep, R_1);
+ return r_minimised(ref1, ref2, items, scalep, R_1_IGNORE);
+}
+
+
+double stat_r1_zero(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep)
+{
+ return r_minimised(ref1, ref2, items, scalep, R_1_ZERO);
}
@@ -298,14 +416,107 @@ double stat_r2(const double *ref1, const double *ref2,
}
-double stat_rdiff(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep)
+double stat_rint(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep)
+{
+ return r_minimised(ref1, ref2, items, scalep, R_INT);
+}
+
+
+double stat_rdiff_zero(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep)
+{
+ return r_minimised(ref1, ref2, items, scalep, R_DIFF_ZERO);
+}
+
+
+double stat_rdiff_ignore(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep)
+{
+ return r_minimised(ref1, ref2, items, scalep, R_DIFF_IGNORE);
+}
+
+
+double stat_pearson_i(const double *ref1, const double *ref2,
+ ReflItemList *items)
+{
+ double *vec1, *vec2;
+ int i = 0;
+ int ni = num_items(items);
+ double val;
+
+ vec1 = malloc(ni*sizeof(double));
+ vec2 = malloc(ni*sizeof(double));
+
+ for ( i=0; i<ni; i++ ) {
+
+ double i1, i2;
+ struct refl_item *it;
+ signed int h, k, l;
+
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
+
+ i1 = lookup_intensity(ref1, h, k, l);
+ i2 = lookup_intensity(ref2, h, k, l);
+
+ vec1[i] = i1;
+ vec2[i] = i2;
+
+ printf("%f %f\n", i1, i2);
+
+ }
+
+ val = gsl_stats_correlation(vec1, 1, vec2, 1, i);
+ free(vec1);
+ free(vec2);
+
+ return val;
+}
+
+
+double stat_pearson_f_ignore(const double *ref1, const double *ref2,
+ ReflItemList *items)
{
- return r_minimised(ref1, ref2, items, scalep, R_DIFF);
+ double *vec1, *vec2;
+ int i = 0;
+ int ni = num_items(items);
+ double val;
+
+ vec1 = malloc(ni*sizeof(double));
+ vec2 = malloc(ni*sizeof(double));
+
+ for ( i=0; i<ni; i++ ) {
+
+ double i1, i2, f1, f2;
+ struct refl_item *it;
+ signed int h, k, l;
+
+ it = get_item(items, i);
+ h = it->h; k = it->k; l = it->l;
+
+ i1 = lookup_intensity(ref1, h, k, l);
+ if ( i1 < 0.0 ) continue;
+ f1 = sqrt(i1);
+ i2 = lookup_intensity(ref2, h, k, l);
+ if ( i2 < 0.0 ) continue;
+ f2 = sqrt(i2);
+
+ vec1[i] = f1;
+ vec2[i] = f2;
+
+ }
+
+ val = gsl_stats_correlation(vec1, 1, vec2, 1, i);
+ free(vec1);
+ free(vec2);
+
+ return val;
}
-double stat_pearson(const double *ref1, const double *ref2, ReflItemList *items)
+double stat_pearson_f_zero(const double *ref1, const double *ref2,
+ ReflItemList *items)
{
double *vec1, *vec2;
int i = 0;
diff --git a/src/statistics.h b/src/statistics.h
index 4a795090..a1b99ff7 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -22,20 +22,31 @@
extern double stat_scale_intensity(const double *ref1, const double *ref2,
ReflItemList *items);
-extern double stat_r1(const double *ref1, const double *ref2,
- ReflItemList *items, double *scalep);
+extern double stat_r1_zero(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep);
+extern double stat_r1_ignore(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep);
extern double stat_r2(const double *ref1, const double *ref2,
ReflItemList *items, double *scalep);
-extern double stat_rdiff(const double *ref1, const double *ref2,
+extern double stat_rint(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep);
+
+extern double stat_rdiff_zero(const double *ref1, const double *ref2,
+ ReflItemList *items, double *scalep);
+extern double stat_rdiff_ignore(const double *ref1, const double *ref2,
ReflItemList *items, double *scalep);
extern double stat_riso(const double *ref1, const double *ref2,
ReflItemList *items, double *scalep);
-extern double stat_pearson(const double *ref1, const double *ref2,
- ReflItemList *items);
+extern double stat_pearson_i(const double *ref1, const double *ref2,
+ ReflItemList *items);
+extern double stat_pearson_f_zero(const double *ref1, const double *ref2,
+ ReflItemList *items);
+extern double stat_pearson_f_ignore(const double *ref1, const double *ref2,
+ ReflItemList *items);
#endif /* STATISTICS_H */
generated by cgit v1.2.3 (git 2.25.1) at 2024-07-27 06:10:03 +0000