aboutsummaryrefslogtreecommitdiff
path: root/src/statistics.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/statistics.c')
-rw-r--r--src/statistics.c257
1 files changed, 234 insertions, 23 deletions
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;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-10 06:26:18 +0000