aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@bitwiz.org.uk>2010-05-04 07:38:31 -0700
committerThomas White <taw@bitwiz.org.uk>2010-05-04 08:56:17 -0700
commite447c2efb128823c93358c51a7d8f23636740f68 (patch)
tree68fb800a56d2e06cc6f3610b58a26ecc9bb6cb57 /src
parent07d406ce05760f633f72c2c3683dd3466423321a (diff)
process_hkl: Implement --scale option
Diffstat (limited to 'src')
-rw-r--r--src/compare_hkl.c7
-rw-r--r--src/likelihood.c44
-rw-r--r--src/likelihood.h29
-rw-r--r--src/list_tmp.h2
-rw-r--r--src/process_hkl.c34
-rw-r--r--src/statistics.c43
-rw-r--r--src/statistics.h8
7 files changed, 129 insertions, 38 deletions
diff --git a/src/compare_hkl.c b/src/compare_hkl.c
index bb07d17e..e7752078 100644
--- a/src/compare_hkl.c
+++ b/src/compare_hkl.c
@@ -53,7 +53,6 @@ int main(int argc, char *argv[])
double scale, R;
unsigned int *c1;
unsigned int *c2;
- unsigned int *cjoint;
int i;
int nc1, nc2, ncom;
@@ -134,18 +133,16 @@ int main(int argc, char *argv[])
}
}
- cjoint = new_list_count();
nc1 = 0;
nc2 = 0;
ncom = 0;
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
- cjoint[i] = c1[i] && c2[i];
nc1 += c1[i];
nc2 += c2[i];
- ncom += cjoint[i];
+ ncom += c1[i] && c2[i];
}
STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom);
- R = stat_r2(ref1, ref2, cjoint, IDIM*IDIM*IDIM, &scale);
+ R = stat_r2(ref1, c1, ref2, c2, &scale);
STATUS("R2 = %5.4f %% (scale=%5.2f)\n", R*100.0, scale);
if ( outfile != NULL ) {
diff --git a/src/likelihood.c b/src/likelihood.c
new file mode 100644
index 00000000..e8c62010
--- /dev/null
+++ b/src/likelihood.c
@@ -0,0 +1,44 @@
+/*
+ * likelihood.c
+ *
+ * Likelihood maximisation
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#include "statistics.h"
+#include "utils.h"
+
+void detwin_intensities(const double *model, double *new_pattern,
+ const unsigned int *model_counts,
+ unsigned int *new_counts)
+{
+ /* Placeholder... */
+}
+
+void scale_intensities(const double *model, double *new_pattern,
+ const unsigned int *model_counts,
+ unsigned int *new_counts)
+{
+ double s;
+ unsigned int i;
+
+ s = stat_scale_intensity(model, model_counts, new_pattern, new_counts);
+ printf("%f\n", s);
+
+ /* NaN -> abort */
+ if ( isnan(s) ) return;
+
+ /* Multiply the new pattern up by "s" */
+ for ( i=0; i<LIST_SIZE; i++ ) {
+ new_counts[i] *= s;
+ }
+}
diff --git a/src/likelihood.h b/src/likelihood.h
new file mode 100644
index 00000000..9a773829
--- /dev/null
+++ b/src/likelihood.h
@@ -0,0 +1,29 @@
+/*
+ * likelihood.h
+ *
+ * Likelihood maximisation
+ *
+ * (c) 2006-2010 Thomas White <taw@physics.org>
+ *
+ * Part of CrystFEL - crystallography with a FEL
+ *
+ */
+
+#ifdef HAVE_CONFIG_H
+#include <config.h>
+#endif
+
+#ifndef LIKELIHOOD_H
+#define LIKELIHOOD_H
+
+
+extern void detwin_intensities(const double *model, double *new_pattern,
+ const unsigned int *model_counts,
+ unsigned int *new_counts);
+
+extern void scale_intensities(const double *model, double *new_pattern,
+ const unsigned int *model_counts,
+ unsigned int *new_counts);
+
+
+#endif /* LIKELIHOOD_H */
diff --git a/src/list_tmp.h b/src/list_tmp.h
index 2b8f703d..749d7e7a 100644
--- a/src/list_tmp.h
+++ b/src/list_tmp.h
@@ -18,6 +18,8 @@
*
*/
+#include <stdio.h>
+
#define ERROR_T(...) fprintf(stderr, __VA_ARGS__)
static inline void LABEL(integrate)(TYPE *ref, signed int h,
diff --git a/src/process_hkl.c b/src/process_hkl.c
index dfb84ee1..d09fca22 100644
--- a/src/process_hkl.c
+++ b/src/process_hkl.c
@@ -66,6 +66,8 @@ static void show_help(const char *s)
" --detwin Correlate each new pattern with the current\n"
" model and choose the best fitting out of the\n"
" allowable twins.\n"
+" --scale Scale each pattern for best fit with the current\n"
+" model.\n"
);
}
@@ -141,8 +143,9 @@ static void write_RvsQ(const char *name, double *ref, double *trueref,
}
-static void process_reflections(double *ref, double *trueref,
- unsigned int *counts, unsigned int n_patterns,
+static void process_reflections(double *ref, unsigned int *counts,
+ double *trueref, unsigned int *truecounts,
+ unsigned int n_patterns,
UnitCell *cell, int do_rvsq, int do_zoneaxis)
{
int j;
@@ -158,7 +161,7 @@ static void process_reflections(double *ref, double *trueref,
}
mean_counts = (double)ctot/nmeas;
- R = stat_r2(ref, trueref, counts, LIST_SIZE, &scale);
+ R = stat_r2(ref, counts, trueref, truecounts, &scale);
STATUS("%8u: R=%5.2f%% (sf=%7.4e) mean meas/refl=%5.2f,"
" %i reflections measured\n",
n_patterns, R*100.0, scale, mean_counts, nmeas);
@@ -236,10 +239,12 @@ int main(int argc, char *argv[])
int config_zoneaxis = 0;
int config_sum = 0;
int config_detwin = 0;
+ int config_scale = 0;
char *intfile = NULL;
double *new_pattern = NULL;
unsigned int *new_counts = NULL;
unsigned int n_total_patterns;
+ unsigned int *truecounts = NULL;
/* Long options */
const struct option longopts[] = {
@@ -254,6 +259,7 @@ int main(int argc, char *argv[])
{"compare-with", 0, NULL, 'c'},
{"sum", 0, &config_sum, 1},
{"detwin", 0, &config_detwin, 1},
+ {"scale", 0, &config_scale, 1},
{0, 0, NULL, 0}
};
@@ -313,8 +319,9 @@ int main(int argc, char *argv[])
}
if ( intfile != NULL ) {
+ truecounts = new_list_count();
STATUS("Comparing against '%s'\n", intfile);
- trueref = read_reflections(intfile, NULL);
+ trueref = read_reflections(intfile, truecounts);
free(intfile);
} else {
trueref = NULL;
@@ -369,19 +376,27 @@ int main(int argc, char *argv[])
continue;
}
+ /* Detwin before scaling */
if ( config_detwin ) {
detwin_intensities(model, new_pattern,
model_counts, new_counts);
}
+ /* Scale if requested */
+ if ( config_scale ) {
+ scale_intensities(model, new_pattern,
+ model_counts, new_counts);
+ }
+
/* Start of second or later pattern */
merge_pattern(model, new_pattern, model_counts,
new_counts, config_maxonly, config_sum);
if (config_every && (n_patterns % config_every == 0)) {
- process_reflections(model, trueref,
- model_counts, n_patterns,
- cell, config_rvsq,
+ process_reflections(model, model_counts,
+ trueref, truecounts,
+ n_patterns, cell,
+ config_rvsq,
config_zoneaxis);
}
@@ -412,8 +427,9 @@ int main(int argc, char *argv[])
fclose(fh);
if ( trueref != NULL ) {
- process_reflections(model, trueref, model_counts, n_patterns,
- cell, config_rvsq, config_zoneaxis);
+ process_reflections(model, model_counts, trueref, truecounts,
+ n_patterns, cell, config_rvsq,
+ config_zoneaxis);
}
if ( output != NULL ) {
diff --git a/src/statistics.c b/src/statistics.c
index be499814..406f1f56 100644
--- a/src/statistics.c
+++ b/src/statistics.c
@@ -18,26 +18,26 @@
#include <stdlib.h>
#include "statistics.h"
+#include "utils.h"
-/* By what (best-fitted) factor must the list "list" be multiplied by,
- * if it were to be merged with "target"? */
-static double stat_scale_intensity(double *obs, double *calc, unsigned int *c,
- int size)
+double stat_scale_intensity(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2)
{
double top = 0.0;
double bot = 0.0;
int i;
/* Start from i=1 -> skip central beam */
- for ( i=1; i<size; i++ ) {
+ for ( i=1; i<LIST_SIZE; i++ ) {
- if ( c[i] > 0 ) {
- double obsi;
- obsi = obs[i] / (double)c[i];
- top += obsi * calc[i];
- bot += calc[i] * calc[i];
- } /* else reflection not measured so don't worry about it */
+ if ( c1[i] && c2[i] ) {
+ double i1, i2;
+ i1 = ref1[i] / (double)c1[i];
+ i2 = ref2[i] / (double)c2[i];
+ top += i1 * i2;
+ bot += i2 * i2;
+ } /* else reflection not common so don't worry about it */
}
@@ -45,28 +45,27 @@ static double stat_scale_intensity(double *obs, double *calc, unsigned int *c,
}
-double stat_r2(double *obs, double *calc, unsigned int *c, int size,
- double *scalep)
+double stat_r2(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2, double *scalep)
{
double top = 0.0;
double bot = 0.0;
double scale;
int i;
- scale = stat_scale_intensity(obs, calc, c, size);
+ scale = stat_scale_intensity(ref1, c1, ref2, c2);
*scalep = scale;
/* Start from i=1 -> skip central beam */
- for ( i=1; i<size; i++ ) {
-
- if ( c[i] > 0 ) {
+ for ( i=1; i<LIST_SIZE; i++ ) {
- double obsi;
+ if ( c1[i] && c2[i] ) {
- obsi = obs[i] / (double)c[i];
- obsi = obsi / scale;
+ double i1, i2;
+ i1 = ref1[i] / (scale*(double)c1[i]);
+ i2 = ref2[i] / (scale*(double)c2[i]);
- top += pow(fabs(obsi - calc[i]), 2.0);
- bot += pow(obsi, 2.0);
+ top += pow(fabs(i1 - i2), 2.0);
+ bot += pow(i1, 2.0);
} /* else reflection not measured so don't worry about it */
diff --git a/src/statistics.h b/src/statistics.h
index 5db7b045..fe7f57e1 100644
--- a/src/statistics.h
+++ b/src/statistics.h
@@ -16,7 +16,11 @@
#ifndef STATISTICS_H
#define STATISTICS_H
-double stat_r2(double *obs, double *calc, unsigned int *c, int size,
- double *scalep);
+extern double stat_scale_intensity(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2);
+
+extern double stat_r2(const double *ref1, const unsigned int *c1,
+ const double *ref2, const unsigned int *c2,
+ double *scalep);
#endif /* STATISTICS_H */