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

process_hkl: Implement --scale option

parent 07d406ce
Hide whitespace changes
Inline Side-by-side
src/compare_hkl.c
View file @ e447c2ef
......@@ -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 ) {
......
src/likelihood.c 0 → 100644
View file @ e447c2ef
/*
* 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;
}
}
src/likelihood.h 0 → 100644
View file @ e447c2ef
/*
* 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 */
src/list_tmp.h
View file @ e447c2ef
......@@ -18,6 +18,8 @@
*
*/
#include <stdio.h>
#define ERROR_T(...) fprintf(stderr, __VA_ARGS__)
static inline void LABEL(integrate)(TYPE *ref, signed int h,
......
src/process_hkl.c
View file @ e447c2ef
......@@ -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 ) {
......
src/statistics.c
View file @ e447c2ef
......@@ -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 */
......
src/statistics.h
View file @ e447c2ef
......@@ -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 */
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