Commit 6a542235 authored by Thomas White's avatar Thomas White
Browse files

Reduce the scope of "count"

Lists of counts had pervaded every corner of CrystFEL, being used as markers
for the presence of reflections.  Now we have a better way of doing this,
the ReflItemList, and few parts of the suite apart from process_hkl have any
business knowing how many observations were made of a particular reflection.
parent 6a476e01
......@@ -49,13 +49,9 @@ int main(int argc, char *argv[])
char *outfile = NULL;
char *afile = NULL;
char *bfile = NULL;
signed int h, k, l;
double scale, R2, Rmerge, pearson;
unsigned int *c1;
unsigned int *c2;
int i;
int nc1, nc2, ncom;
unsigned int *outcounts;
int i, ncom;
ReflItemList *i1, *i2, *icommon;
/* Long options */
const struct option longopts[] = {
......@@ -94,72 +90,53 @@ int main(int argc, char *argv[])
bfile = strdup(argv[optind]);
cell = load_cell_from_pdb("molecule.pdb");
c1 = new_list_count();
ref1 = read_reflections(afile, c1, NULL);
ref1 = new_list_intensity();
i1 = read_reflections(afile, ref1, NULL, NULL);
if ( ref1 == NULL ) {
ERROR("Couldn't open file '%s'\n", afile);
return 1;
}
c2 = new_list_count();
ref2 = read_reflections(bfile, c2, NULL);
ref2 = new_list_intensity();
i2 = read_reflections(bfile, ref2, NULL, NULL);
if ( ref2 == NULL ) {
ERROR("Couldn't open file '%s'\n", bfile);
return 1;
}
out = new_list_intensity();
outcounts = new_list_count();
/* Knock out the zero beam, in case it's present */
set_count(c1, 0, 0, 0, 0);
set_count(c2, 0, 0, 0, 0);
/* Find common reflections */
icommon = intersection_items(i1, i2);
ncom = num_items(icommon);
/* Divide by number of counts, since we're not interested in them */
divide_down(ref1, c1);
divide_down(ref2, c2);
/* List for output scale factor map */
out = new_list_intensity();
for ( h=-INDMAX; h<INDMAX; h++ ) {
for ( k=-INDMAX; k<INDMAX; k++ ) {
for ( l=-INDMAX; l<INDMAX; l++ ) {
for ( i=0; i<ncom; i++ ) {
double i1, i2;
unsigned int c1s, c2s;
struct refl_item *it;
signed int h, k, l;
c1s = lookup_count(c1, h, k, l);
c2s = lookup_count(c2, h, k, l);
it = get_item(icommon, i);
h = it->h; k = it->k; l = it->l;
i1 = lookup_intensity(ref1, h, k, l);
i2 = lookup_intensity(ref2, h, k, l);
if ( c1s && c2s ) {
if ( (i1 != 0.0) && (i2 != 0.0) ) {
set_intensity(out, h, k, l,
(i1/(double)c1s)/i2/(double)c2s);
set_count(outcounts, h, k, l, 1);
}
}
set_intensity(out, h, k, l, i1/i2);
}
}
}
nc1 = 0;
nc2 = 0;
ncom = 0;
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
nc1 += c1[i];
nc2 += c2[i];
ncom += c1[i] && c2[i];
}
STATUS("%i,%i reflections: %i in common\n", nc1, nc2, ncom);
R2 = stat_r2(ref1, c1, ref2, c2, &scale);
STATUS("%i,%i reflections: %i in common\n",
num_items(i1), num_items(i2), ncom);
R2 = stat_r2(ref1, ref2, icommon, &scale);
STATUS("R2 = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale);
Rmerge = stat_rmerge(ref1, c1, ref2, c2, &scale);
Rmerge = stat_rmerge(ref1, ref2, icommon, &scale);
STATUS("Rmerge = %5.4f %% (scale=%5.2e)\n", Rmerge*100.0, scale);
pearson = stat_pearson(ref1, c1, ref2, c2);
pearson = stat_pearson(ref1, ref2, icommon);
STATUS("Pearson r = %5.4f\n", pearson);
if ( outfile != NULL ) {
write_reflections(outfile, outcounts, out, NULL, 0, cell, 1);
write_reflections(outfile, icommon, out, NULL, NULL, cell);
}
return 0;
......
......@@ -342,8 +342,7 @@ void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
/* Setup the OpenCL stuff, create buffers, load the structure factor table */
struct gpu_context *setup_gpu(int no_sfac, struct image *image,
const double *intensities,
const unsigned int *counts)
const double *intensities)
{
struct gpu_context *gctx;
cl_uint nplat;
......@@ -404,11 +403,7 @@ struct gpu_context *setup_gpu(int no_sfac, struct image *image,
intensities_ptr = malloc(intensities_size);
if ( intensities != NULL ) {
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
if ( counts[i] == 1 ) {
intensities_ptr[i] = intensities[i];
} else {
intensities_ptr[i] = 1.0e20;
}
intensities_ptr[i] = intensities[i];
}
} else {
for ( i=0; i<IDIM*IDIM*IDIM; i++ ) {
......
......@@ -26,8 +26,7 @@ struct gpu_context;
extern void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
int na, int nb, int nc, UnitCell *ucell);
extern struct gpu_context *setup_gpu(int no_sfac, struct image *image,
const double *intensities,
const unsigned int *counts);
const double *intensities);
extern void cleanup_gpu(struct gpu_context *gctx);
#else
......@@ -40,8 +39,7 @@ static void get_diffraction_gpu(struct gpu_context *gctx, struct image *image,
}
static struct gpu_context *setup_gpu(int no_sfac, struct image *image,
const double *intensities,
const unsigned int *counts)
const double *intensities)
{
return NULL;
}
......
......@@ -71,13 +71,11 @@ static double lattice_factor(struct rvec q, double ax, double ay, double az,
static double interpolate_linear(const double *ref,
const unsigned int *counts,
float hd, signed int k, signed int l)
{
signed int h;
double val1, val2;
float f;
unsigned int c1, c2;
h = (signed int)hd;
if ( hd < 0.0 ) h -= 1;
......@@ -86,30 +84,15 @@ static double interpolate_linear(const double *ref,
val1 = lookup_intensity(ref, h, k, l);
val2 = lookup_intensity(ref, h+1, k, l);
c1 = lookup_count(counts, h, k, l);
c2 = lookup_count(counts, h+1, k, l);
if ( c1 == 0 ) {
ERROR("Needed intensity for %i %i %i, but don't have it.\n",
h, k, l);
return 1.0e20;
}
if ( c2 == 0 ) {
ERROR("Needed intensity for %i %i %i, but don't have it.\n",
h+1, k, l);
return 1.0e20;
}
val1 = val1 / (double)c1;
val2 = val2 / (double)c2;
val1 = val1;
val2 = val2;
return (1.0-f)*val1 + f*val2;
}
static double interpolate_bilinear(const double *ref,
const unsigned int *counts,
float hd, float kd, signed int l)
{
signed int k;
......@@ -121,15 +104,14 @@ static double interpolate_bilinear(const double *ref,
f = kd - (float)k;
assert(f >= 0.0);
val1 = interpolate_linear(ref, counts, hd, k, l);
val2 = interpolate_linear(ref, counts, hd, k+1, l);
val1 = interpolate_linear(ref, hd, k, l);
val2 = interpolate_linear(ref, hd, k+1, l);
return (1.0-f)*val1 + f*val2;
}
static double interpolate_intensity(const double *ref,
const unsigned int *counts,
float hd, float kd, float ld)
{
signed int l;
......@@ -141,15 +123,14 @@ static double interpolate_intensity(const double *ref,
f = ld - (float)l;
assert(f >= 0.0);
val1 = interpolate_bilinear(ref, counts, hd, kd, l);
val2 = interpolate_bilinear(ref, counts, hd, kd, l+1);
val1 = interpolate_bilinear(ref, hd, kd, l);
val2 = interpolate_bilinear(ref, hd, kd, l+1);
return (1.0-f)*val1 + f*val2;
}
static double complex interpolate_phased_linear(const double *ref,
const unsigned int *counts,
const double *phases,
float hd,
signed int k, signed int l)
......@@ -157,7 +138,6 @@ static double complex interpolate_phased_linear(const double *ref,
signed int h;
double val1, val2;
float f;
unsigned int c1, c2;
double ph1, ph2;
double re1, re2, im1, im2;
double re, im;
......@@ -169,25 +149,11 @@ static double complex interpolate_phased_linear(const double *ref,
val1 = lookup_intensity(ref, h, k, l);
val2 = lookup_intensity(ref, h+1, k, l);
c1 = lookup_count(counts, h, k, l);
c2 = lookup_count(counts, h+1, k, l);
ph1 = lookup_phase(phases, h, k, l);
ph2 = lookup_phase(phases, h+1, k, l);
if ( c1 == 0 ) {
ERROR("Needed intensity for %i %i %i, but don't have it.\n",
h, k, l);
return 1.0e20;
}
if ( c2 == 0 ) {
ERROR("Needed intensity for %i %i %i, but don't have it.\n",
h+1, k, l);
return 1.0e20;
}
val1 = val1 / (double)c1;
val2 = val2 / (double)c2;
val1 = val1;
val2 = val2;
/* Calculate real and imaginary parts */
re1 = val1 * cos(ph1);
......@@ -203,7 +169,6 @@ static double complex interpolate_phased_linear(const double *ref,
static double complex interpolate_phased_bilinear(const double *ref,
const unsigned int *counts,
const double *phases,
float hd, float kd,
signed int l)
......@@ -217,15 +182,14 @@ static double complex interpolate_phased_bilinear(const double *ref,
f = kd - (float)k;
assert(f >= 0.0);
val1 = interpolate_phased_linear(ref, counts, phases, hd, k, l);
val2 = interpolate_phased_linear(ref, counts, phases, hd, k+1, l);
val1 = interpolate_phased_linear(ref, phases, hd, k, l);
val2 = interpolate_phased_linear(ref, phases, hd, k+1, l);
return (1.0-f)*val1 + f*val2;
}
static double interpolate_phased_intensity(const double *ref,
const unsigned int *counts,
const double *phases,
float hd, float kd, float ld)
{
......@@ -238,16 +202,15 @@ static double interpolate_phased_intensity(const double *ref,
f = ld - (float)l;
assert(f >= 0.0);
val1 = interpolate_phased_bilinear(ref, counts, phases, hd, kd, l);
val2 = interpolate_phased_bilinear(ref, counts, phases, hd, kd, l+1);
val1 = interpolate_phased_bilinear(ref, phases, hd, kd, l);
val2 = interpolate_phased_bilinear(ref, phases, hd, kd, l+1);
return cabs((1.0-f)*val1 + f*val2);
}
/* Look up the structure factor for the nearest Bragg condition */
static double molecule_factor(const double *intensities,
const unsigned int *counts, const double *phases,
static double molecule_factor(const double *intensities,const double *phases,
struct rvec q,
double ax, double ay, double az,
double bx, double by, double bz,
......@@ -267,18 +230,13 @@ static double molecule_factor(const double *intensities,
h = (signed int)rint(hd);
k = (signed int)rint(kd);
l = (signed int)rint(ld);
if ( lookup_count(counts, h, k, l) == 0 ) {
ERROR("Needed intensity for %i %i %i, but don't have it.\n",
h, k, l);
return 1.0e20;
}
r = lookup_intensity(intensities, h, k, l);
break;
case GRADIENT_INTERPOLATE :
r = interpolate_intensity(intensities, counts, hd, kd, ld);
r = interpolate_intensity(intensities, hd, kd, ld);
break;
case GRADIENT_PHASED :
r = interpolate_phased_intensity(intensities, counts, phases,
r = interpolate_phased_intensity(intensities, phases,
hd, kd, ld);
break;
default:
......@@ -291,7 +249,7 @@ static double molecule_factor(const double *intensities,
double water_diffraction(struct rvec q, double en,
double beam_r, double water_r)
double beam_r, double water_r)
{
double complex fH, fO;
double s, modq;
......@@ -390,9 +348,8 @@ double get_tt(struct image *image, unsigned int xs, unsigned int ys)
void get_diffraction(struct image *image, int na, int nb, int nc,
const double *intensities, const unsigned int *counts,
const double *phases, UnitCell *cell, int do_water,
GradientMethod m)
const double *intensities, const double *phases,
UnitCell *cell, int do_water, GradientMethod m)
{
unsigned int xs, ys;
double ax, ay, az;
......@@ -447,7 +404,7 @@ void get_diffraction(struct image *image, int na, int nb, int nc,
I_molecule = 1.0e10;
} else {
I_molecule = molecule_factor(intensities,
counts, phases, q,
phases, q,
ax,ay,az,
bx,by,bz,cx,cy,cz,
m);
......
......@@ -26,8 +26,7 @@ typedef enum {
} GradientMethod;
extern void get_diffraction(struct image *image, int na, int nb, int nc,
const double *intensities,
const unsigned int *counts, const double *phases,
const double *intensities,const double *phases,
UnitCell *cell, int do_water, GradientMethod m);
extern struct rvec get_q(struct image *image, unsigned int xs, unsigned int ys,
unsigned int sampling, float *ttp, float k);
......
......@@ -38,8 +38,6 @@ static void show_help(const char *s)
" --poisson Simulate Poisson samples.\n"
" --twin Generate twinned data.\n"
" -o, --output=<filename> Output filename (default: stdout).\n"
" --zone-axis Generate hk0 intensities only (and add\n"
" Synth2D-style header.\n"
" -i, --intensities=<file> Read intensities from file instead of\n"
" calculating them from scratch. You might use\n"
" this if you need to apply noise or twinning.\n"
......@@ -48,37 +46,6 @@ static void show_help(const char *s)
}
static int template_reflections(const char *filename, unsigned int *counts)
{
char *rval;
FILE *fh;
fh = fopen(filename, "r");
if ( fh == NULL ) {
return 1;
}
do {
char line[1024];
int r;
signed int h, k, l;
rval = fgets(line, 1023, fh);
r = sscanf(line, "%i %i %i", &h, &k, &l);
if ( r != 3 ) continue;
set_count(counts, h, k, l, 1);
} while ( rval != NULL );
fclose(fh);
return 0;
}
/* Apply Poisson noise to all reflections */
static void noisify_reflections(double *ref)
{
......@@ -111,13 +78,12 @@ int main(int argc, char *argv[])
char *template = NULL;
int config_noisify = 0;
int config_twin = 0;
int config_za = 0;
char *output = NULL;
unsigned int *counts;
unsigned int *cts;
char *input = NULL;
signed int h, k, l;
char *filename = NULL;
ReflItemList *input_items;
ReflItemList *write_items;
/* Long options */
const struct option longopts[] = {
......@@ -126,7 +92,6 @@ int main(int argc, char *argv[])
{"poisson", 0, &config_noisify, 1},
{"output", 1, NULL, 'o'},
{"twin", 0, &config_twin, 1},
{"zone-axis", 0, &config_za, 1},
{"intensities", 1, NULL, 'i'},
{"pdb", 1, NULL, 'p'},
{0, 0, NULL, 0}
......@@ -170,41 +135,18 @@ int main(int argc, char *argv[])
}
mol = load_molecule(filename);
cts = new_list_count();
phases = new_list_intensity(); /* "intensity" type used for phases */
phases = new_list_phase();
if ( input == NULL ) {
input_items = new_items();
ideal_ref = get_reflections(mol, eV_to_J(1790.0), 1/(0.05e-9),
cts, phases);
phases, input_items);
} else {
ideal_ref = read_reflections(input, cts, phases);
ideal_ref = new_list_intensity();
phases = new_list_phase();
input_items = read_reflections(input, ideal_ref, phases, NULL);
free(input);
}
counts = new_list_count();
if ( template != NULL ) {
if ( template_reflections(template, counts) != 0 ) {
ERROR("Failed to template reflections.\n");
return 1;
}
} else {
/* No template? Then only mark reflections which were
* calculated. */
for ( h=-INDMAX; h<=INDMAX; h++ ) {
for ( k=-INDMAX; k<=INDMAX; k++ ) {
for ( l=-INDMAX; l<=INDMAX; l++ ) {
unsigned int c;
c = lookup_count(cts, h, k, l);
set_count(counts, h, k, l, c);
}
}
}
}
if ( config_noisify ) noisify_reflections(ideal_ref);
if ( config_twin ) {
......@@ -239,8 +181,24 @@ int main(int argc, char *argv[])
}
write_reflections(output, counts, ideal_ref, phases,
config_za, mol->cell, 1);
if ( template ) {
/* Write out only reflections which are in the template
* (and which we have in the input) */
ReflItemList *template_items;
template_items = read_reflections(template, NULL, NULL, NULL);
write_items = intersection_items(input_items, template_items);
delete_items(template_items);
} else {
/* Write out all reflections */
write_items = new_items();
union_items(write_items, input_items);
}
write_reflections(output, write_items, ideal_ref, phases, NULL,
mol->cell);
delete_items(input_items);
delete_items(write_items);
return 0;
}
......@@ -65,7 +65,6 @@ struct process_args
struct detector *det;
IndexingMethod indm;
const double *intensities;
const unsigned int *counts;
struct gpu_context *gctx;
/* Thread control and output */
......@@ -225,19 +224,18 @@ static struct image *get_simage(struct image *template, int alternate)
static void simulate_and_write(struct image *simage, struct gpu_context **gctx,
const double *intensities,
const unsigned int *counts, UnitCell *cell)
const double *intensities, UnitCell *cell)
{
/* Set up GPU if necessary */
if ( (gctx != NULL) && (*gctx == NULL) ) {
*gctx = setup_gpu(0, simage, intensities, counts);
*gctx = setup_gpu(0, simage, intensities);
}
if ( (gctx != NULL) && (*gctx != NULL) ) {
get_diffraction_gpu(*gctx, simage, 24, 24, 40, cell);
} else {
get_diffraction(simage, 24, 24, 40,
intensities, counts, NULL, cell, 0,
intensities, NULL, cell, 0,
GRADIENT_MOSAIC);
}
record_image(simage, 0);
......@@ -270,7 +268,6 @@ static struct process_result process_image(struct process_args *pargs)
int config_polar = pargs->config_polar;
IndexingMethod indm = pargs->indm;
const double *intensities = pargs->intensities;
const unsigned int *counts = pargs->counts;
struct gpu_context *gctx = pargs->gctx;
image.features = NULL;
......@@ -368,11 +365,11 @@ static struct process_result process_image(struct process_args *pargs)
if ( config_gpu ) {
pthread_mutex_lock(pargs->gpu_mutex);
simulate_and_write(simage, &gctx, intensities,
counts, image.indexed_cell);
image.indexed_cell);
pthread_mutex_unlock(pargs->gpu_mutex);
} else {
simulate_and_write(simage, NULL, intensities,
counts, image.indexed_cell);
image.indexed_cell);
}
}
......@@ -469,7 +466,6 @@ int main(int argc, char *argv[])
UnitCell *cell;
double *intensities = NULL;
char *intfile = NULL;
unsigned int *counts;
char *pdb = NULL;
char *prefix = NULL;
int nthreads = 1;
......@@ -568,11 +564,11 @@ int main(int argc, char *argv[])
free(filename);
if ( intfile != NULL ) {
counts = new_list_count();
intensities = read_reflections(intfile, counts, NULL);
ReflItemList *items;
items = read_reflections(intfile, intensities, NULL, NULL);
delete_items(items);
} else {
intensities = NULL;
counts = NULL;
}
if ( pdb == NULL ) {
......@@ -671,7 +667,6 @@ int main(int argc, char *argv[])
pargs->det = det;
pargs->indm = indm;
pargs->intensities = intensities;
pargs->counts = counts;
pargs->gctx = gctx;
pargs->id = i;
pthread_mutex_lock(&pargs->control_mutex);
......
......@@ -18,17 +18,17 @@
#include "utils.h"
void scale_intensities(const double *model, double *new_pattern,
const unsigned int *model_counts,
ReflItemList *items, double f0, int f0_valid)
void scale_intensities(const double *model, ReflItemList *model_items,
double *new_pattern, ReflItemList *new_items,
double f0, int f0_valid)
{