From 6264ff6e0ac13fcc22afcb4ce9f98a920444a92c Mon Sep 17 00:00:00 2001
From: Thomas White <taw@physics.org>
Date: Mon, 12 Jul 2010 10:06:47 +0200
Subject: compare_hkl: Add Pearson correlation

---
 src/compare_hkl.c |  4 +++-
 src/statistics.c  | 36 ++++++++++++++++++++++++++++++++++++
 src/statistics.h  |  3 +++
 3 files changed, 42 insertions(+), 1 deletion(-)

(limited to 'src')

diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index 3e313733..0dffd3a0 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -50,7 +50,7 @@ int main(int argc, char *argv[])
 	char *afile = NULL;
 	char *bfile = NULL;
 	signed int h, k, l;
-	double scale, R2, Rmerge;
+	double scale, R2, Rmerge, pearson;
 	unsigned int *c1;
 	unsigned int *c2;
 	int i;
@@ -145,6 +145,8 @@ int main(int argc, char *argv[])
 	STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale);
 	Rmerge = stat_rmerge(ref1, c1, ref2, c2, &scale);
 	STATUS("Rmerge = %5.4f %% (scale=%5.2e)\n", Rmerge*100.0, scale);
+	pearson = stat_pearson(ref1, c1, ref2, c2);
+	STATUS("Pearson r = %5.4f\n", pearson);
 
 	if ( outfile != NULL ) {
 		write_reflections(outfile, NULL, out, NULL, 1, cell, 1);
diff --git a/src/statistics.c b/src/statistics.c
index da3d95c4..f3fb21c0 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -18,6 +18,7 @@
 #include <stdlib.h>
 #include <gsl/gsl_errno.h>
 #include <gsl/gsl_min.h>
+#include <gsl/gsl_statistics.h>
 
 #include "statistics.h"
 #include "utils.h"
@@ -249,3 +250,38 @@ double stat_r2(const double *ref1, const unsigned int *c1,
 {
 	return r_minimised(ref1, c1, ref2, c2, scalep, R_2);
 }
+
+
+double stat_pearson(const double *ref1, const unsigned int *c1,
+                    const double *ref2, const unsigned int *c2)
+{
+	double vec1[4096];
+	double vec2[4096];
+	signed int h, k, l;
+	int i = 0;
+
+	for ( l=-INDMAX; l<INDMAX; l++ ) {
+	for ( k=-INDMAX; k<INDMAX; k++ ) {
+	for ( h=-INDMAX; h<INDMAX; h++ ) {
+
+		double i1, i2;
+		unsigned int c1s, c2s;
+
+		c1s = lookup_count(c1, h, k, l);
+		c2s = lookup_count(c2, h, k, l);
+
+		i1 = lookup_intensity(ref1, h, k, l);
+		i2 = lookup_intensity(ref2, h, k, l);
+
+		if ( c1s && c2s && (i1>0.0) && (i2>0.0) ) {
+			vec1[i] = sqrt(i1 / (double)c1s);
+			vec2[i] = sqrt(i2 / (double)c2s);
+			i++;
+		}
+
+	}
+	}
+	}
+
+	return gsl_stats_correlation(vec1, 1, vec2, 1, i);
+}
diff --git a/src/statistics.h b/src/statistics.h
index e19d4b29..5e3ca52c 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -27,4 +27,7 @@ extern double stat_rmerge(const double *ref1, const unsigned int *c1,
                           const double *ref2, const unsigned int *c2,
                           double *scalep);
 
+extern double stat_pearson(const double *ref1, const unsigned int *c1,
+                           const double *ref2, const unsigned int *c2);
+
 #endif	/* STATISTICS_H */
-- 
cgit v1.2.3