From 6a5422356c15962726df2261aa53354b0ff12662 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 14 Jul 2010 17:58:55 +0200 Subject: 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. --- src/compare_hkl.c | 69 ++++++------------ src/diffraction-gpu.c | 9 +-- src/diffraction-gpu.h | 6 +- src/diffraction.c | 81 +++++---------------- src/diffraction.h | 3 +- src/get_hkl.c | 94 +++++++----------------- src/indexamajig.c | 21 ++---- src/likelihood.c | 16 ++-- src/likelihood.h | 7 +- src/pattern_sim.c | 15 ++-- src/process_hkl.c | 80 ++++++++++++-------- src/reflections.c | 83 +++++++++------------ src/reflections.h | 14 ++-- src/render_hkl.c | 10 ++- src/sfac.c | 4 +- src/sfac.h | 3 +- src/statistics.c | 197 ++++++++++++++++++++++++++------------------------ src/statistics.h | 22 +++--- src/utils.c | 42 ++++++++--- src/utils.h | 2 + 20 files changed, 345 insertions(+), 433 deletions(-) (limited to 'src') diff --git a/src/compare_hkl.c b/src/compare_hkl.c index e6caee6f..e69f3753 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -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; hh; 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= 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); diff --git a/src/diffraction.h b/src/diffraction.h index c7afaa3e..ffe3a19d 100644 --- a/src/diffraction.h +++ b/src/diffraction.h @@ -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); diff --git a/src/get_hkl.c b/src/get_hkl.c index c98508ec..78d79f76 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -38,8 +38,6 @@ static void show_help(const char *s) " --poisson Simulate Poisson samples.\n" " --twin Generate twinned data.\n" " -o, --output= Output filename (default: stdout).\n" -" --zone-axis Generate hk0 intensities only (and add\n" -" Synth2D-style header.\n" " -i, --intensities= 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; } diff --git a/src/indexamajig.c b/src/indexamajig.c index 3a3dfb56..5df08ad3 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -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); diff --git a/src/likelihood.c b/src/likelihood.c index b6a37994..cb5595f8 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -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) { double s; unsigned int i; - unsigned int *new_counts; + ReflItemList *items; - new_counts = items_to_counts(items); - - s = stat_scale_intensity(model, model_counts, new_pattern, new_counts); + items = intersection_items(model_items, new_items); + s = stat_scale_intensity(model, new_pattern, items); + delete_items(items); if ( f0_valid ) printf("%f %f\n", s, f0); /* NaN -> abort */ @@ -38,6 +38,4 @@ void scale_intensities(const double *model, double *new_pattern, for ( i=0; iop; @@ -134,8 +136,9 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins, } - fom = stat_pearson(trial_ints, trial_counts, - model, model_counts); + intersection = intersection_items(observed, items); + fom = stat_pearson(trial_ints, model, intersection); + delete_items(intersection); free(trial_ints); free(trial_counts); @@ -153,9 +156,9 @@ static int resolve_twin(ReflItemList *items, ReflItemList *twins, } -static void merge_pattern(double *model, const double *new, - unsigned int *model_counts, - ReflItemList *items, int mo, int sum, +static void merge_pattern(double *model, ReflItemList *observed, + const double *new, ReflItemList *items, + unsigned int *model_counts, int mo, ReflItemList *twins, const char *holo, const char *mero) { @@ -163,8 +166,8 @@ static void merge_pattern(double *model, const double *new, int twin; if ( twins != NULL ) { - twin = resolve_twin(items, twins, new, model, - model_counts, holo, mero); + twin = resolve_twin(model, observed, new, items, + twins, holo, mero); } else { twin = 0; } @@ -187,25 +190,23 @@ static void merge_pattern(double *model, const double *new, intensity = lookup_intensity(new, h, k, l); + /* User asked for max only? */ if ( !mo ) { integrate_intensity(model, h, k, l, intensity); - if ( sum ) { - set_count(model_counts, h, k, l, 1); - } else { - integrate_count(model_counts, h, k, l, 1); - } } else { if ( intensity > lookup_intensity(model, h, k, l) ) { set_intensity(model, h, k, l, intensity); } - set_count(model_counts, h, k, l, 1); } + integrate_count(model_counts, h, k, l, 1); + } } -static void merge_all(FILE *fh, double *model, unsigned int *model_counts, +static void merge_all(FILE *fh, double **pmodel, ReflItemList **pobserved, + unsigned int **pcounts, int config_maxonly, int config_scale, int config_sum, int config_stopafter, ReflItemList *twins, const char *holo, const char *mero, @@ -218,6 +219,10 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts, int n_patterns = 0; double *new_pattern = new_list_intensity(); ReflItemList *items = new_items(); + ReflItemList *observed = new_items(); + double *model = new_list_intensity(); + unsigned int *counts = new_list_count(); + int i; do { @@ -244,19 +249,20 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts, /* Scale if requested */ if ( config_scale ) { - scale_intensities(model, new_pattern, - model_counts, items, f0, - f0_valid); + scale_intensities(model, observed, + new_pattern, items, + f0, f0_valid); } /* Start of second or later pattern */ - merge_pattern(model, new_pattern, model_counts, - items, config_maxonly, config_sum, twins, - holo, mero); + merge_pattern(model, observed, new_pattern, items, + counts, config_maxonly, + twins, holo, mero); if ( n_patterns == config_stopafter ) break; n_patterns++; + union_items(observed, items); clear_items(items); progress_bar(n_patterns, n_total_patterns, "Merging"); @@ -292,6 +298,19 @@ static void merge_all(FILE *fh, double *model, unsigned int *model_counts, delete_items(items); free(new_pattern); + if ( !config_sum ) { + for ( i=0; i 0 ) { + model[i] /= (double)counts[i]; + counts[i] = 1; + } + } + } + + *pmodel = model; + *pcounts = counts; + *pobserved = observed; + STATUS("%i patterns had no f0 valid value.\n", n_nof0); } @@ -322,7 +341,7 @@ int main(int argc, char *argv[]) char *output = NULL; FILE *fh; double *model; - unsigned int *model_counts; + unsigned int *counts; UnitCell *cell; int config_maxonly = 0; int config_stopafter = 0; @@ -332,6 +351,7 @@ int main(int argc, char *argv[]) char *sym = NULL; char *pdb = NULL; ReflItemList *twins; + ReflItemList *observed; int i; /* Long options */ @@ -400,9 +420,6 @@ int main(int argc, char *argv[]) sym = strdup("1"); } - model = new_list_intensity(); - model_counts = new_list_count(); - cell = load_cell_from_pdb(pdb); free(pdb); @@ -443,7 +460,7 @@ int main(int argc, char *argv[]) STATUS("There are %i patterns to process\n", n_total_patterns); rewind(fh); - merge_all(fh, model, model_counts, + merge_all(fh, &model, &observed, &counts, config_maxonly, config_scale, config_sum, config_stopafter, twins, holo, sym, n_total_patterns); rewind(fh); @@ -451,13 +468,12 @@ int main(int argc, char *argv[]) fclose(fh); if ( output != NULL ) { - write_reflections(output, model_counts, model, NULL, - 0, cell, 1); + write_reflections(output, observed, model, NULL, counts, cell); } free(sym); free(model); - free(model_counts); + free(counts); free(output); free(cell); diff --git a/src/reflections.c b/src/reflections.c index bd9b446d..0ab7a91e 100644 --- a/src/reflections.c +++ b/src/reflections.c @@ -21,12 +21,12 @@ #include "reflections.h" -void write_reflections(const char *filename, unsigned int *counts, - double *ref, double *phases, int zone_axis, - UnitCell *cell, unsigned int min_counts) +void write_reflections(const char *filename, ReflItemList *items, + double *intensities, double *phases, + unsigned int *counts, UnitCell *cell) { FILE *fh; - signed int h, k, l; + int i; if ( filename == NULL ) { fh = stdout; @@ -40,23 +40,16 @@ void write_reflections(const char *filename, unsigned int *counts, } /* Write spacings and angle if zone axis pattern */ - if ( zone_axis ) { - double a, b, c, alpha, beta, gamma; - cell_get_parameters(cell, &a, &b, &c, &alpha, &beta, &gamma); - fprintf(fh, "a %5.3f nm\n", - (0.5/resolution(cell, 1, 0, 0))*1e9); - fprintf(fh, "b %5.3f nm\n", - (0.5/resolution(cell, 0, 1, 0))*1e9); - fprintf(fh, "angle %5.3f deg\n", rad2deg(alpha)); - fprintf(fh, "scale 10\n"); - } else { - fprintf(fh, " h k l I phase sigma(I) " - " 1/d(nm^-1) counts\n"); - } + fprintf(fh, " h k l I phase sigma(I) " + " 1/d(nm^-1) counts\n"); + + for ( i=0; ih; k = it->k; l = it->l; int N; double intensity, s; @@ -64,14 +57,12 @@ void write_reflections(const char *filename, unsigned int *counts, if ( counts ) { N = lookup_count(counts, h, k, l); - if ( N < min_counts ) continue; } else { N = 1; } - if ( zone_axis && (l != 0) ) continue; - /* Divide measured intensity by the number of counts */ - intensity = lookup_intensity(ref, h, k, l) / N; + intensity = lookup_intensity(intensities, h, k, l); + if ( phases != NULL ) { double p; p = lookup_phase(phases, h, k, l); @@ -90,19 +81,26 @@ void write_reflections(const char *filename, unsigned int *counts, fprintf(fh, "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i\n", h, k, l, intensity, ph, 0.0, s/1.0e9, N); - } - } } fclose(fh); } -double *read_reflections(const char *filename, unsigned int *counts, - double *phases) +/* Read reflections from file. Returns the list of reflections successfully + * read in. "intensities", "phases" and "counts" are lists which will be + * populated with the values read from the file. Existing values in either list + * will be overwritten if the reflection is read from the file, but other values + * will be left intact. + * + * "intensities", "phases" or "counts" can be NULL, if you don't need them. + */ +ReflItemList *read_reflections(const char *filename, + double *intensities, double *phases, + unsigned int *counts) { - double *ref; FILE *fh; char *rval; + ReflItemList *items; fh = fopen(filename, "r"); if ( fh == NULL ) { @@ -110,7 +108,7 @@ double *read_reflections(const char *filename, unsigned int *counts, return NULL; } - ref = new_list_intensity(); + items = new_items(); do { @@ -147,24 +145,24 @@ double *read_reflections(const char *filename, unsigned int *counts, continue; } - set_intensity(ref, h, k, l, intensity); + add_item(items, h, k, l); + + if ( intensities != NULL ) { + set_intensity(intensities, h, k, l, intensity); + } if ( phases != NULL ) { ph = atof(phs); set_phase(phases, h, k, l, ph); } if ( counts != NULL ) { set_count(counts, h, k, l, cts); - /* In this case, the intensity must be multiplied up - * because other parts of the program will try to - * divide it down. */ - set_intensity(ref, h, k, l, intensity*(double)cts); } } while ( rval != NULL ); fclose(fh); - return ref; + return items; } @@ -188,16 +186,3 @@ double *ideal_intensities(double complex *sfac) return ref; } - - -void divide_down(double *intensities, unsigned int *counts) -{ - int i; - - for ( i=0; i 0 ) { - intensities[i] /= (double)counts[i]; - counts[i] = 1; - } - } -} diff --git a/src/reflections.h b/src/reflections.h index 84946f1a..de5feb72 100644 --- a/src/reflections.h +++ b/src/reflections.h @@ -18,18 +18,18 @@ #include "cell.h" +#include "utils.h" -extern void write_reflections(const char *filename, unsigned int *counts, - double *ref, double *phases, int zone_axis, - UnitCell *cell, unsigned int min_counts); +extern void write_reflections(const char *filename, ReflItemList *items, + double *intensities, double *phases, + unsigned int *counts, UnitCell *cell); -extern double *read_reflections(const char *filename, unsigned int *counts, - double *phases); +extern ReflItemList *read_reflections(const char *filename, + double *intensities, double *phases, + unsigned int *counts); extern double *ideal_intensities(double complex *sfac); -extern void divide_down(double *intensities, unsigned int *counts); - #endif /* REFLECTIONS_H */ diff --git a/src/render_hkl.c b/src/render_hkl.c index 60153652..829ceb4e 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -128,10 +128,10 @@ static void render_za(UnitCell *cell, double *ref, unsigned int *c, switch ( wght ) { case WGHT_I : - val = lookup_intensity(ref, h, k, 0) / (float)ct; + val = lookup_intensity(ref, h, k, 0); break; case WGHT_SQRTI : - val = lookup_intensity(ref, h, k, 0) / (float)ct; + val = lookup_intensity(ref, h, k, 0); val = (val>0.0) ? sqrt(val) : 0.0; break; case WGHT_COUNTS : @@ -290,7 +290,6 @@ int main(int argc, char *argv[]) UnitCell *cell; char *infile; double *ref; - unsigned int *cts; int config_povray = 0; int config_zoneaxis = 0; int config_sqrt = 0; @@ -303,6 +302,7 @@ int main(int argc, char *argv[]) int wght; int colscale; char *cscale = NULL; + unsigned int *cts; /* Long options */ const struct option longopts[] = { @@ -417,8 +417,10 @@ int main(int argc, char *argv[]) ERROR("Couldn't load unit cell from %s\n", pdb); return 1; } + ref = new_list_intensity(); cts = new_list_count(); - ref = read_reflections(infile, cts, NULL); + ReflItemList *items = read_reflections(infile, ref, NULL, cts); + delete_items(items); if ( ref == NULL ) { ERROR("Couldn't open file '%s'\n", infile); return 1; diff --git a/src/sfac.c b/src/sfac.c index a8b94051..b5ae7ace 100644 --- a/src/sfac.c +++ b/src/sfac.c @@ -491,7 +491,7 @@ void free_molecule(struct molecule *mol) double *get_reflections(struct molecule *mol, double en, double res, - unsigned int *counts, double *phases) + double *phases, ReflItemList *items) { double *reflections; double asx, asy, asz; @@ -556,7 +556,7 @@ double *get_reflections(struct molecule *mol, double en, double res, set_intensity(reflections, h, k, l, pow(cabs(F), 2.0)); if ( phases != NULL ) set_phase(phases, h, k, l, carg(F)); - set_count(counts, h, k, l, 1); + if ( items != NULL ) add_item(items, h, k, l); } progress_bar((k+INDMAX)+IDIM*(h+INDMAX), IDIM*IDIM-1, diff --git a/src/sfac.h b/src/sfac.h index ea30c83e..115dcb5f 100644 --- a/src/sfac.h +++ b/src/sfac.h @@ -19,6 +19,7 @@ #include #include "cell.h" +#include "utils.h" #define MAX_ATOMS (128*1024) @@ -61,7 +62,7 @@ extern double complex get_sfac(const char *n, double s, double en); extern struct molecule *load_molecule(const char *filename); extern void free_molecule(struct molecule *mol); extern double *get_reflections(struct molecule *mol, double en, double res, - unsigned int *counts, double *phases); + double *phases, ReflItemList *items); #endif /* SFAC_H */ diff --git a/src/statistics.c b/src/statistics.c index 889ade51..1acbb718 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -26,10 +26,9 @@ struct r_params { const double *ref1; - const unsigned int *c1; const double *ref2; - const unsigned int *c2; - int fom; + ReflItemList *items; /* Which reflections to use */ + int fom; /* Which FoM to use (see the enum just below) */ }; enum { @@ -38,23 +37,31 @@ enum { }; -double stat_scale_intensity(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2) +/* Return the least squares optimal scaling factor when comparing intensities. + * ref1,ref2 are the two intensity lists to compare. "items" is a ReflItemList + * containing the reflections which should be taken into account. + */ +double stat_scale_intensity(const double *ref1, const double *ref2, + ReflItemList *items) { double top = 0.0; double bot = 0.0; int i; - /* Start from i=1 -> skip central beam */ - for ( i=1; ih; k = it->k; l = it->l; + + i1 = lookup_intensity(ref1, h, k, l); + i2 = lookup_intensity(ref2, h, k, l); - 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 */ + top += i1 * i2; + bot += i2 * i2; } @@ -62,29 +69,36 @@ double stat_scale_intensity(const double *ref1, const unsigned int *c1, } -double stat_scale_sqrti(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2) +/* Return the least squares optimal scaling factor when comparing the square + * roots of the intensities (i.e. one approximation to the structure factor + * moduli). + * ref1,ref2 are the two intensity lists to compare (they contain intensities, + * not square rooted intensities). "items" is a ReflItemList containing the + * reflections which should be taken into account. + */ +static double stat_scale_sqrti(const double *ref1, const double *ref2, + ReflItemList *items) { double top = 0.0; double bot = 0.0; int i; - /* Start from i=1 -> skip central beam */ - for ( i=1; ih; k = it->k; l = it->l; - f1 = sqrt(ref1[i]) / (double)c1[i]; - f2 = sqrt(ref2[i]) / (double)c2[i]; - - top += f1 * f2; - bot += f2 * f2; + i1 = lookup_intensity(ref1, h, k, l); + f1 = i1 > 0.0 ? sqrt(i1) : 0.0; + i2 = lookup_intensity(ref2, h, k, l); + f2 = i2 > 0.0 ? sqrt(i2) : 0.0; - } /* else reflection not common so don't worry about it */ + top += i1 * i2; + bot += i2 * i2; } @@ -92,27 +106,27 @@ double stat_scale_sqrti(const double *ref1, const unsigned int *c1, } -static double internal_r2(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2, - double scale) +static double internal_r2(const double *ref1, const double *ref2, + ReflItemList *items, double scale) { double top = 0.0; double bot = 0.0; int i; - /* Start from i=1 -> skip central beam */ - for ( i=1; ih; k = it->k; l = it->l; - top += pow(i1 - i2, 2.0); - bot += pow(i1, 2.0); + i1 = scale * lookup_intensity(ref1, h, k, l); + i2 = lookup_intensity(ref2, h, k, l); - } /* else reflection not measured so don't worry about it */ + top += pow(i1 - i2, 2.0); + bot += pow(i1, 2.0); } @@ -120,30 +134,29 @@ static double internal_r2(const double *ref1, const unsigned int *c1, } -static double internal_rmerge(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2, - double scale) +static double internal_rmerge(const double *ref1, const double *ref2, + ReflItemList *items, double scale) { double top = 0.0; double bot = 0.0; int i; - /* Start from i=1 -> skip central beam */ - for ( i=1; ih; k = it->k; l = it->l; - f1 = sqrt(ref1[i]) / (scale*(double)c1[i]); - f2 = sqrt(ref2[i]) / (double)c2[i]; - - top += fabs(f1 - f2); - bot += f1 + f2; + i1 = lookup_intensity(ref1, h, k, l); + f1 = i1 > 0.0 ? sqrt(i1) : 0.0; + i2 = lookup_intensity(ref2, h, k, l); + f2 = i2 > 0.0 ? sqrt(i2) : 0.0; - } /* else reflection not measured so don't worry about it */ + top += fabs(f1 - f2); + bot += f1 + f2; } @@ -157,11 +170,9 @@ static double calc_r(double scale, void *params) switch ( rp->fom) { case R_MERGE : - return internal_rmerge(rp->ref1, rp->c1, - rp->ref2, rp->c2, scale); + return internal_rmerge(rp->ref1, rp->ref2, rp->items, scale); case R_2 : - return internal_r2(rp->ref1, rp->c1, - rp->ref2, rp->c2, scale); + return internal_r2(rp->ref1, rp->ref2, rp->items, scale); } ERROR("No such FoM!\n"); @@ -169,9 +180,8 @@ static double calc_r(double scale, void *params) } -static double r_minimised(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2, - double *scalep, int fom) +static double r_minimised(const double *ref1, const double *ref2, + ReflItemList *items, double *scalep, int fom) { gsl_function F; gsl_min_fminimizer *s; @@ -182,8 +192,7 @@ static double r_minimised(const double *ref1, const unsigned int *c1, rp.ref1 = ref1; rp.ref2 = ref2; - rp.c1 = c1; - rp.c2 = c2; + rp.items = items; rp.fom = fom; F.function = &calc_r; @@ -194,10 +203,10 @@ static double r_minimised(const double *ref1, const unsigned int *c1, /* Initial guess */ switch ( fom ) { case R_MERGE : - scale = stat_scale_sqrti(ref1, c1, ref2, c2); + scale = stat_scale_sqrti(ref1, ref2, items); break; case R_2 : - scale = stat_scale_intensity(ref1, c1, ref2, c2); + scale = stat_scale_intensity(ref1, ref2, items); break; } //STATUS("Initial scale factor estimate: %5.2e\n", scale); @@ -236,52 +245,52 @@ static double r_minimised(const double *ref1, const unsigned int *c1, } -double stat_rmerge(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2, - double *scalep) +double stat_rmerge(const double *ref1, const double *ref2, + ReflItemList *items, double *scalep) { - return r_minimised(ref1, c1, ref2, c2, scalep, R_MERGE); + return r_minimised(ref1, ref2, items, scalep, R_MERGE); } -double stat_r2(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2, - double *scalep) +double stat_r2(const double *ref1, const double *ref2, + ReflItemList *items, double *scalep) { - return r_minimised(ref1, c1, ref2, c2, scalep, R_2); + return r_minimised(ref1, ref2, items, scalep, R_2); } -double stat_pearson(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2) +double stat_pearson(const double *ref1, const double *ref2, ReflItemList *items) { - double vec1[4096]; - double vec2[4096]; - signed int h, k, l; + double *vec1, *vec2; int i = 0; + int ni = num_items(items); + double val; - for ( l=-INDMAX; lh; k = it->k; l = it->l; i1 = lookup_intensity(ref1, h, k, l); + f1 = i1 > 0.0 ? sqrt(i1) : 0.0; i2 = lookup_intensity(ref2, h, k, l); + f2 = i2 > 0.0 ? sqrt(i2) : 0.0; - if ( c1s && c2s && (i1>0.0) && (i2>0.0) ) { - vec1[i] = sqrt(i1 / (double)c1s); - vec2[i] = sqrt(i2 / (double)c2s); - i++; - } + vec1[i] = f1; + vec2[i] = f2; } - } - } - return gsl_stats_correlation(vec1, 1, vec2, 1, i); + val = gsl_stats_correlation(vec1, 1, vec2, 1, i); + free(vec1); + free(vec2); + + return val; } diff --git a/src/statistics.h b/src/statistics.h index 5e3ca52c..c6c2f09c 100644 --- a/src/statistics.h +++ b/src/statistics.h @@ -16,18 +16,20 @@ #ifndef STATISTICS_H #define STATISTICS_H -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); +#include "utils.h" -extern double stat_rmerge(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2, - double *scalep); +extern double stat_scale_intensity(const double *ref1, const double *ref2, + ReflItemList *items); + +extern double stat_rmerge(const double *ref1, const double *ref2, + ReflItemList *items, double *scalep); + +extern double stat_r2(const double *ref1, const double *ref2, + ReflItemList *items, double *scalep); + +extern double stat_pearson(const double *ref1, const double *ref2, + ReflItemList *items); -extern double stat_pearson(const double *ref1, const unsigned int *c1, - const double *ref2, const unsigned int *c2); #endif /* STATISTICS_H */ diff --git a/src/utils.c b/src/utils.c index 21ba882f..0f20435e 100644 --- a/src/utils.c +++ b/src/utils.c @@ -387,24 +387,23 @@ int num_items(const ReflItemList *items) } -unsigned int *items_to_counts(ReflItemList *items) +void union_op_items(ReflItemList *items, ReflItemList *newi) { - unsigned int *c; - int i; + int n, i; - c = new_list_count(); + n = num_items(newi); + for ( i=0; ih, r->k, r->l, 1); - } + struct refl_item *r = get_item(newi, i); + if ( find_op(items, r->op) ) continue; + + add_item_with_op(items, r->h, r->k, r->l, r->op); - return c; + } } -void union_op_items(ReflItemList *items, ReflItemList *newi) +void union_items(ReflItemList *items, ReflItemList *newi) { int n, i; @@ -412,9 +411,28 @@ void union_op_items(ReflItemList *items, ReflItemList *newi) for ( i=0; iop) ) continue; + if ( find_item(items, r->h, r->k, r->l) ) continue; add_item_with_op(items, r->h, r->k, r->l, r->op); } } + + +ReflItemList *intersection_items(ReflItemList *i1, ReflItemList *i2) +{ + int n, i; + ReflItemList *res = new_items(); + + n = num_items(i1); + for ( i=0; ih, r->k, r->l) ) { + add_item_with_op(res, r->h, r->k, r->l, r->op); + } + + } + + return res; +} diff --git a/src/utils.h b/src/utils.h index a7433d28..f7e98469 100644 --- a/src/utils.h +++ b/src/utils.h @@ -200,6 +200,8 @@ extern struct refl_item *get_item(ReflItemList *items, int i); extern int num_items(const ReflItemList *items); extern unsigned int *items_to_counts(ReflItemList *items); extern void union_op_items(ReflItemList *items, ReflItemList *newi); +extern void union_items(ReflItemList *items, ReflItemList *newi); +extern ReflItemList *intersection_items(ReflItemList *i1, ReflItemList *i2); /* ------------------------------ Message macros ---------------------------- */ -- cgit v1.2.3