From b8ac0faea1abb2e3033140dcd157160bfa98d251 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 15 Mar 2011 15:56:26 +0100 Subject: WIP on new stream format in process_hkl --- Makefile.am | 14 +- Makefile.in | 43 ++-- src/process_hkl.c | 570 ++++++++++++---------------------------------------- src/reflist-utils.c | 72 +++++++ src/reflist-utils.h | 30 +++ src/reflist.c | 43 ++++ src/reflist.h | 6 + src/stream.c | 128 ++++++------ src/stream.h | 4 + 9 files changed, 389 insertions(+), 521 deletions(-) create mode 100644 src/reflist-utils.c create mode 100644 src/reflist-utils.h diff --git a/Makefile.am b/Makefile.am index ce2d06d7..f7437acc 100644 --- a/Makefile.am +++ b/Makefile.am @@ -42,10 +42,10 @@ tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c src/utils.c \ endif src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \ - src/cell.c src/utils.c src/reflections.c \ + src/cell.c src/utils.c \ src/symmetry.c src/stream.c src/beam-parameters.c \ src/thread-pool.c src/image.c src/detector.c \ - src/hdf5-file.c src/reflist.c + src/hdf5-file.c src/reflist.c src/reflist-utils.c src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \ src/cell.c src/image.c src/peaks.c src/index.c \ @@ -53,7 +53,8 @@ src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c src/utils.c \ src/sfac.c src/dirax.c src/mosflm.c \ src/reflections.c src/symmetry.c \ src/geometry.c src/thread-pool.c \ - src/beam-parameters.c src/reflist.c src/stream.c + src/beam-parameters.c src/reflist.c src/stream.c \ + src/reflist-utils.c if HAVE_OPENCL src_indexamajig_SOURCES += src/diffraction-gpu.c src/cl-utils.c endif @@ -99,18 +100,19 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \ src/geometry.c src/reflections.c src/stream.c \ src/thread-pool.c src/beam-parameters.c \ src/symmetry.c src/post-refinement.c \ - src/hrs-scaling.c src/reflist.c + src/hrs-scaling.c src/reflist.c src/reflist-utils.c if BUILD_CUBEIT src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \ src/detector.c src/render.c src/filters.c src/image.c \ - src/symmetry.c src/stream.c src/thread-pool.c src/reflist.c + src/symmetry.c src/stream.c src/thread-pool.c \ + src/reflist.c src/reflist-utils.c endif src_estimate_background_SOURCES = src/estimate_background.c src/stream.c \ src/utils.c src/cell.c src/thread-pool.c \ src/image.c src/detector.c src/hdf5-file.c \ - src/reflist.c + src/reflist.c src/reflist-utils.c tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \ src/utils.c diff --git a/Makefile.in b/Makefile.in index b65b60a3..402094df 100644 --- a/Makefile.in +++ b/Makefile.in @@ -107,7 +107,7 @@ src_compare_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am__src_cubeit_SOURCES_DIST = src/cubeit.c src/cell.c src/hdf5-file.c \ src/utils.c src/detector.c src/render.c src/filters.c \ src/image.c src/symmetry.c src/stream.c src/thread-pool.c \ - src/reflist.c + src/reflist.c src/reflist-utils.c @BUILD_CUBEIT_TRUE@am_src_cubeit_OBJECTS = src/cubeit.$(OBJEXT) \ @BUILD_CUBEIT_TRUE@ src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) \ @BUILD_CUBEIT_TRUE@ src/utils.$(OBJEXT) src/detector.$(OBJEXT) \ @@ -115,7 +115,8 @@ am__src_cubeit_SOURCES_DIST = src/cubeit.c src/cell.c src/hdf5-file.c \ @BUILD_CUBEIT_TRUE@ src/image.$(OBJEXT) src/symmetry.$(OBJEXT) \ @BUILD_CUBEIT_TRUE@ src/stream.$(OBJEXT) \ @BUILD_CUBEIT_TRUE@ src/thread-pool.$(OBJEXT) \ -@BUILD_CUBEIT_TRUE@ src/reflist.$(OBJEXT) +@BUILD_CUBEIT_TRUE@ src/reflist.$(OBJEXT) \ +@BUILD_CUBEIT_TRUE@ src/reflist-utils.$(OBJEXT) src_cubeit_OBJECTS = $(am_src_cubeit_OBJECTS) src_cubeit_LDADD = $(LDADD) src_cubeit_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -124,7 +125,7 @@ am_src_estimate_background_OBJECTS = \ src/utils.$(OBJEXT) src/cell.$(OBJEXT) \ src/thread-pool.$(OBJEXT) src/image.$(OBJEXT) \ src/detector.$(OBJEXT) src/hdf5-file.$(OBJEXT) \ - src/reflist.$(OBJEXT) + src/reflist.$(OBJEXT) src/reflist-utils.$(OBJEXT) src_estimate_background_OBJECTS = \ $(am_src_estimate_background_OBJECTS) src_estimate_background_LDADD = $(LDADD) @@ -154,8 +155,8 @@ am__src_indexamajig_SOURCES_DIST = src/indexamajig.c src/hdf5-file.c \ src/filters.c src/diffraction.c src/detector.c src/sfac.c \ src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \ src/geometry.c src/thread-pool.c src/beam-parameters.c \ - src/reflist.c src/stream.c src/diffraction-gpu.c \ - src/cl-utils.c + src/reflist.c src/stream.c src/reflist-utils.c \ + src/diffraction-gpu.c src/cl-utils.c @HAVE_OPENCL_TRUE@am__objects_1 = src/diffraction-gpu.$(OBJEXT) \ @HAVE_OPENCL_TRUE@ src/cl-utils.$(OBJEXT) am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \ @@ -166,7 +167,8 @@ am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \ src/mosflm.$(OBJEXT) src/reflections.$(OBJEXT) \ src/symmetry.$(OBJEXT) src/geometry.$(OBJEXT) \ src/thread-pool.$(OBJEXT) src/beam-parameters.$(OBJEXT) \ - src/reflist.$(OBJEXT) src/stream.$(OBJEXT) $(am__objects_1) + src/reflist.$(OBJEXT) src/stream.$(OBJEXT) \ + src/reflist-utils.$(OBJEXT) $(am__objects_1) src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS) src_indexamajig_LDADD = $(LDADD) src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -177,7 +179,7 @@ am_src_partialator_OBJECTS = src/partialator.$(OBJEXT) \ src/stream.$(OBJEXT) src/thread-pool.$(OBJEXT) \ src/beam-parameters.$(OBJEXT) src/symmetry.$(OBJEXT) \ src/post-refinement.$(OBJEXT) src/hrs-scaling.$(OBJEXT) \ - src/reflist.$(OBJEXT) + src/reflist.$(OBJEXT) src/reflist-utils.$(OBJEXT) src_partialator_OBJECTS = $(am_src_partialator_OBJECTS) src_partialator_LDADD = $(LDADD) src_partialator_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -205,11 +207,11 @@ src_powder_plot_LDADD = $(LDADD) src_powder_plot_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_src_process_hkl_OBJECTS = src/process_hkl.$(OBJEXT) \ src/sfac.$(OBJEXT) src/statistics.$(OBJEXT) src/cell.$(OBJEXT) \ - src/utils.$(OBJEXT) src/reflections.$(OBJEXT) \ - src/symmetry.$(OBJEXT) src/stream.$(OBJEXT) \ - src/beam-parameters.$(OBJEXT) src/thread-pool.$(OBJEXT) \ - src/image.$(OBJEXT) src/detector.$(OBJEXT) \ - src/hdf5-file.$(OBJEXT) src/reflist.$(OBJEXT) + src/utils.$(OBJEXT) src/symmetry.$(OBJEXT) \ + src/stream.$(OBJEXT) src/beam-parameters.$(OBJEXT) \ + src/thread-pool.$(OBJEXT) src/image.$(OBJEXT) \ + src/detector.$(OBJEXT) src/hdf5-file.$(OBJEXT) \ + src/reflist.$(OBJEXT) src/reflist-utils.$(OBJEXT) src_process_hkl_OBJECTS = $(am_src_process_hkl_OBJECTS) src_process_hkl_LDADD = $(LDADD) src_process_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -612,17 +614,17 @@ src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c \ @HAVE_OPENCL_TRUE@ src/image.c src_process_hkl_SOURCES = src/process_hkl.c src/sfac.c src/statistics.c \ - src/cell.c src/utils.c src/reflections.c \ + src/cell.c src/utils.c \ src/symmetry.c src/stream.c src/beam-parameters.c \ src/thread-pool.c src/image.c src/detector.c \ - src/hdf5-file.c src/reflist.c + src/hdf5-file.c src/reflist.c src/reflist-utils.c src_indexamajig_SOURCES = src/indexamajig.c src/hdf5-file.c \ src/utils.c src/cell.c src/image.c src/peaks.c src/index.c \ src/filters.c src/diffraction.c src/detector.c src/sfac.c \ src/dirax.c src/mosflm.c src/reflections.c src/symmetry.c \ src/geometry.c src/thread-pool.c src/beam-parameters.c \ - src/reflist.c src/stream.c $(am__append_6) + src/reflist.c src/stream.c src/reflist-utils.c $(am__append_6) @BUILD_HDFSEE_TRUE@src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c \ @BUILD_HDFSEE_TRUE@ src/hdf5-file.c src/utils.c src/image.c src/filters.c \ @BUILD_HDFSEE_TRUE@ src/thread-pool.c src/detector.c @@ -661,16 +663,17 @@ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \ src/geometry.c src/reflections.c src/stream.c \ src/thread-pool.c src/beam-parameters.c \ src/symmetry.c src/post-refinement.c \ - src/hrs-scaling.c src/reflist.c + src/hrs-scaling.c src/reflist.c src/reflist-utils.c @BUILD_CUBEIT_TRUE@src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \ @BUILD_CUBEIT_TRUE@ src/detector.c src/render.c src/filters.c src/image.c \ -@BUILD_CUBEIT_TRUE@ src/symmetry.c src/stream.c src/thread-pool.c src/reflist.c +@BUILD_CUBEIT_TRUE@ src/symmetry.c src/stream.c src/thread-pool.c \ +@BUILD_CUBEIT_TRUE@ src/reflist.c src/reflist-utils.c src_estimate_background_SOURCES = src/estimate_background.c src/stream.c \ src/utils.c src/cell.c src/thread-pool.c \ src/image.c src/detector.c src/hdf5-file.c \ - src/reflist.c + src/reflist.c src/reflist-utils.c tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \ src/utils.c @@ -851,6 +854,8 @@ src/stream.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/reflist.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) +src/reflist-utils.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) src/cubeit$(EXEEXT): $(src_cubeit_OBJECTS) $(src_cubeit_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/cubeit$(EXEEXT) $(AM_V_CCLD)$(LINK) $(src_cubeit_OBJECTS) $(src_cubeit_LDADD) $(LIBS) @@ -977,6 +982,7 @@ mostlyclean-compile: -rm -f src/powder_plot.$(OBJEXT) -rm -f src/process_hkl.$(OBJEXT) -rm -f src/reflections.$(OBJEXT) + -rm -f src/reflist-utils.$(OBJEXT) -rm -f src/reflist.$(OBJEXT) -rm -f src/render.$(OBJEXT) -rm -f src/render_hkl.$(OBJEXT) @@ -1024,6 +1030,7 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/powder_plot.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/process_hkl.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflections.Po@am__quote@ +@AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist-utils.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/reflist.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render.Po@am__quote@ @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/render_hkl.Po@am__quote@ diff --git a/src/process_hkl.c b/src/process_hkl.c index 8df2dae1..da3ee3f8 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -3,7 +3,7 @@ * * Assemble and process FEL Bragg intensities * - * (c) 2006-2010 Thomas White + * (c) 2006-2011 Thomas White * * Part of CrystFEL - crystallography with a FEL * @@ -24,10 +24,12 @@ #include "utils.h" #include "statistics.h" #include "sfac.h" -#include "reflections.h" +#include "reflist-utils.h" #include "symmetry.h" #include "stream.h" #include "beam-parameters.h" +#include "reflist.h" +#include "image.h" /* Number of divisions for intensity histograms */ @@ -43,7 +45,7 @@ static void show_help(const char *s) " -h, --help Display this help message.\n" " -i, --input= Specify input filename (\"-\" for stdin).\n" " -o, --output= Specify output filename for merged intensities\n" -" (don't specify for no output).\n" +" Default: processed.hkl).\n" " -p, --pdb= PDB file to use (default: molecule.pdb).\n" " -b, --beam= Get beam parameters from file (used for sigmas).\n" "\n" @@ -66,13 +68,10 @@ static void show_help(const char *s) " --scale Scale each pattern for best fit with the current\n" " model.\n" " -y, --symmetry= Merge according to point group .\n" -" -a, --input-symmetry= Specify the apparent (input) symmetry.\n" " --reference= Compare against intensities from when\n" " scaling or resolving ambiguities.\n" " The symmetry of the reference list must be the\n" " same as that given with '-y'.\n" -" --outstream= Write an annotated version of the input stream\n" -" to .\n" ); } @@ -119,175 +118,55 @@ static void plot_histogram(double *vals, int n) } -/* Note "holo" needn't actually be a holohedral point group, if you want to try - * something strange like resolving from a low-symmetry group into an even - * lower symmetry one. - */ -static ReflItemList *get_twin_possibilities(const char *holo, const char *mero) -{ - ReflItemList *test_items; - ReflItemList *twins; - int np; - - np = num_general_equivs(holo) / num_general_equivs(mero); - - test_items = new_items(); - - /* Some arbitrarily chosen reflections which can't be special - * reflections in any point group, i.e. lots of odd numbers, - * prime numbers and so on. There's probably an analytical - * way of working these out, but this will do. */ - add_item(test_items, 1, 2, 3); - add_item(test_items, 3, 7, 13); - add_item(test_items, 5, 2, 1); - - twins = get_twins(test_items, holo, mero); - delete_items(test_items); - - /* Idiot check. Wouldn't be necessary if I could prove that the above - * set of arbitrarily chosen reflections were always general. */ - if ( num_items(twins) != np ) { - ERROR("Whoops! Couldn't find all the twinning possiblities.\n"); - abort(); - } - - return twins; -} - - -static int resolve_twin(const double *model, ReflItemList *observed, - const double *patt, ReflItemList *items, - ReflItemList *twins, const char *holo, const char *mero) -{ - int n, i; - double best_fom = 0.0; - int best_op = 0; - - n = num_items(twins); - - for ( i=0; iop; - - for ( j=0; jh, r->k, r->l, &h, &k, &l, - holo, op); - get_asymm(h, k, l, &h, &k, &l, mero); - - set_intensity(trial_ints, h, k, l, - lookup_intensity(patt, r->h, r->k, r->l)); - set_count(trial_counts, h, k, l, 1); - - } - - intersection = intersection_items(observed, items); - fom = stat_pearson_i(trial_ints, model, intersection); - delete_items(intersection); - - free(trial_ints); - free(trial_counts); - - //printf(" %f", fom); - if ( fom > best_fom ) { - best_fom = fom; - best_op = op; - } - - } - //printf("\n"); - - return best_op; -} - - -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, - ReflItemList *reference, const double *reference_i, - double *hist_vals, - signed int hist_h, signed int hist_k, - signed int hist_l, int *hist_n, double *devs, - double *means, FILE *outfh) +static void merge_pattern(RefList *model, RefList *new, int max_only, + const char *sym, + double *hist_vals, signed int hist_h, + signed int hist_k, signed int hist_l, int *hist_n) { - int i; - int twin; - ReflItemList *sym_items = new_items(); - - if ( twins != NULL ) { - if ( reference != NULL ) { - twin = resolve_twin(reference_i, reference, new, items, - twins, holo, mero); - } else { - twin = resolve_twin(model, observed, new, items, - twins, holo, mero); - } - } else { - twin = 0; - } - - if ( outfh != NULL ) { - fprintf(outfh, "twin=%i\n", twin); - } + Reflection *refl; + RefListIterator *iter; - for ( i=0; ih; - ks = item->k; - ls = item->l; - - /* Transform into correct side of the twin law. - * "twin" is always zero if no de-twinning is performed. */ - get_general_equiv(hs, ks, ls, &h, &k, &l, holo, twin); + get_indices(refl, &h, &k, &l); + model_version = find_refl(model, h, k, l); + if ( model_version == NULL ) { + model_version = add_refl(model, h, k, l); + } /* Put into the asymmetric unit for the target group */ - get_asymm(h, k, l, &h, &k, &l, mero); + get_asymm(h, k, l, &h, &k, &l, sym); /* Read the intensity from the original location * (i.e. before screwing around with symmetry) */ - intensity = lookup_intensity(new, hs, ks, ls); + intensity = get_intensity(refl); + + /* Get the current model intensity */ + model_int = get_intensity(model_version); /* User asked for max only? */ - if ( !mo ) { - integrate_intensity(model, h, k, l, intensity); + if ( !max_only ) { + set_int(model_version, model_int + intensity); } else { - if ( intensity > lookup_intensity(model, h, k, l) ) { - set_intensity(model, h, k, l, intensity); + if ( intensity > get_intensity(model_version) ) { + set_int(model_version, intensity); } } - if ( devs != NULL ) { - double m; - m = lookup_intensity(means, h, k, l); - integrate_intensity(devs, h, k, l, - pow(intensity-m, 2.0)); - } - - if ( !find_item(sym_items, h, k, l) ) { - add_item(sym_items, h, k, l); - } + double dev = get_sum_squared_dev(model_version); + set_sum_squared_dev(model_version, + dev + pow(intensity-model_int, 2.0)); - /* Increase count */ - integrate_count(model_counts, h, k, l, 1); + /* Increase redundancy */ + int cur_redundancy = get_redundancy(model_version); + set_redundancy(model_version, cur_redundancy+1); if ( hist_vals != NULL ) { int p = *hist_n; @@ -298,11 +177,6 @@ static void merge_pattern(double *model, ReflItemList *observed, } } - - /* Dump the reflections in this pattern into the overall list */ - union_items(observed, sym_items); - - delete_items(sym_items); } @@ -314,34 +188,36 @@ enum { }; -static void scale_intensities(const double *model, ReflItemList *model_items, - double *new_pattern, ReflItemList *new_items, - double f0, int f0_valid, const char *sym) +static void scale_intensities(RefList *model, RefList *new, const char *sym) { double s; double top = 0.0; double bot = 0.0; - unsigned int i; const int scaling = SCALE_INTPERBRAGG; + Reflection *refl; + RefListIterator *iter; + Reflection *model_version; - for ( i=0; ih, it->k, it->l, sym, - &hu, &ku, &lu); + model_version = find_refl(model, h, k, l); + if ( model_version == NULL ) continue; + + get_asymm(h, k, l, &hu, &ku, &lu, sym); - i1 = lookup_intensity(model, hu, ku, lu); - i2 = lookup_intensity(new_pattern, it->h, it->k, it->l); + i1 = get_intensity(model_version); + i2 = get_intensity(refl); /* Calculate LSQ estimate of scaling factor */ top += i1 * i2; @@ -352,7 +228,7 @@ static void scale_intensities(const double *model, ReflItemList *model_items, case SCALE_CONSTINT : /* Sum up the intensity in the pattern */ - i2 = lookup_intensity(new_pattern, it->h, it->k, it->l); + i2 = get_intensity(refl); top += i2; break; @@ -360,7 +236,7 @@ static void scale_intensities(const double *model, ReflItemList *model_items, case SCALE_INTPERBRAGG : /* Sum up the intensity in the pattern */ - i2 = lookup_intensity(new_pattern, it->h, it->k, it->l); + i2 = get_intensity(refl); top += i2; bot += 1.0; @@ -381,162 +257,108 @@ static void scale_intensities(const double *model, ReflItemList *model_items, s = 1000.0 / (top/bot); break; } - //if ( f0_valid ) printf("%f %f\n", s, f0); /* Multiply the new pattern up by "s" */ - for ( i=0; i0) ) { + /* Get data from next chunk */ + rval = read_chunk(fh, &image); + if ( rval ) continue; - /* Assume a default I0 if we don't have one by now */ - if ( config_scale && !f0_valid ) { - n_nof0++; - f0 = 1.0; - } + n_patterns++; + + if ( (image.reflections != NULL) && (image.indexed_cell) ) { /* Scale if requested */ if ( config_scale ) { - if ( reference == NULL ) { - scale_intensities(model, observed, - new_pattern, items, - f0, f0_valid, mero); - } else { - scale_intensities(reference_i, - reference, - new_pattern, items, - f0, f0_valid, mero); - } + scale_intensities(model, image.reflections, + sym); } - /* Start of second or later pattern */ - merge_pattern(model, observed, new_pattern, items, - counts, config_maxonly, - twins, holo, mero, - reference, reference_i, - hist_vals, hist_h, hist_k, hist_l, - hist_i, devs, means, outfh); - - n_patterns++; - if ( n_patterns == config_stopafter ) { - progress_bar(n_total_patterns, n_total_patterns, - "Merging"); - break; - } - progress_bar(n_patterns, n_total_patterns, "Merging"); + merge_pattern(model, image.reflections, config_maxonly, + sym, hist_vals, hist_h, hist_k, hist_l, + hist_i); - /* Reset for the next pattern */ - clear_items(items); - - f0_valid = 0; + n_used++; } - if ( outfh != NULL ) { - fprintf(outfh, "%s", line); - } + reflist_free(image.reflections); + image_feature_list_free(image.features); + cell_free(image.indexed_cell); - if ( strncmp(line, "f0 = ", 5) == 0 ) { - r = sscanf(line, "f0 = %f", &f0); - if ( r != 1 ) { - f0 = 1.0; - f0_valid = 0; - continue; - } - f0_valid = 1; - } + progress_bar(n_patterns, n_total_patterns-config_startafter, + "Merging"); - r = sscanf(line, "%i %i %i %f", &h, &k, &l, &intensity); - if ( r != 4 ) continue; + } while ( rval == 0 ); - /* Not interested in the central beam */ - if ( (h==0) && (k==0) && (l==0) ) continue; + /* Divide by counts to get mean intensity if necessary */ + if ( !config_sum && !config_maxonly ) { - /* The same raw indices (before mapping into the asymmetric - * unit should not turn up twice in one pattern. */ - if ( find_item(items, h, k, l) != 0 ) { - ERROR("More than one measurement for %i %i %i in" - " pattern number %i\n", h, k, l, n_patterns); - } - set_intensity(new_pattern, h, k, l, intensity); + Reflection *refl; + RefListIterator *iter; - /* NB: This list contains raw indices, before working out - * where they belong in the asymmetric unit. */ - add_item(items, h, k, l); + for ( refl = first_refl(model, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - } while ( rval != NULL ); + double intensity = get_intensity(refl); + double sum_squared_dev = get_sum_squared_dev(refl); + int red = get_redundancy(refl); - delete_items(items); - free(new_pattern); + set_int(refl, intensity / (double)red); + set_sum_squared_dev(refl, sum_squared_dev/(double)red); - /* Calculate mean intensity if necessary */ - if ( !config_sum && !config_maxonly ) { - for ( i=0; i 0 ) { - model[i] /= (double)counts[i]; - } } + } - *pmodel = model; - *pcounts = counts; - *pobserved = observed; + /* Calculate ESDs */ + for ( refl = first_refl(model, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + double sum_squared_dev = get_sum_squared_dev(refl); + int red = get_redundancy(refl); + + set_esd_intensity(refl, sum_squared_dev/(double)red); + + } - STATUS("%i patterns had no f0 valid value.\n", n_nof0); + STATUS("%i of the patterns could be used.\n", n_used); } @@ -546,35 +368,22 @@ int main(int argc, char *argv[]) char *filename = NULL; char *output = NULL; FILE *fh; - double *model; - unsigned int *counts; + RefList *model; UnitCell *cell = NULL; int config_maxonly = 0; int config_startafter = 0; int config_stopafter = 0; int config_sum = 0; int config_scale = 0; - int config_rmerge = 0; unsigned int n_total_patterns; char *sym = NULL; - char *in_sym = NULL; char *pdb = NULL; - ReflItemList *twins; - ReflItemList *observed; - int i; char *histo = NULL; signed int hist_h, hist_k, hist_l; double *hist_vals = NULL; int hist_i; - char *outstream = NULL; - char *reference = NULL; - ReflItemList *reference_items; - double *reference_i; - FILE *outfh = NULL; struct beam_params *beam = NULL; int space_for_hist = 0; - double *devs; - double *esds; /* Long options */ const struct option longopts[] = { @@ -588,18 +397,14 @@ int main(int argc, char *argv[]) {"sum", 0, &config_sum, 1}, {"scale", 0, &config_scale, 1}, {"symmetry", 1, NULL, 'y'}, - {"input-symmetry", 1, NULL, 'a'}, {"pdb", 1, NULL, 'p'}, {"histogram", 1, NULL, 'g'}, - {"rmerge", 0, &config_rmerge, 1}, - {"outstream", 1, NULL, 2}, - {"reference", 1, NULL, 'r'}, {"beam", 1, NULL, 'b'}, {0, 0, NULL, 0} }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:e:r:o:p:y:g:f:a:b:", + while ((c = getopt_long(argc, argv, "hi:e:o:p:y:g:f:b:", longopts, NULL)) != -1) { switch (c) { @@ -635,18 +440,6 @@ int main(int argc, char *argv[]) histo = strdup(optarg); break; - case 'r' : - reference = strdup(optarg); - break; - - case 2 : - outstream = strdup(optarg); - break; - - case 'a' : - in_sym = strdup(optarg); - break; - case 'b' : beam = get_beam_parameters(optarg); if ( beam == NULL ) { @@ -665,17 +458,15 @@ int main(int argc, char *argv[]) } - if ( config_sum && config_rmerge ) { - ERROR("Options --sum and --rmerge do not make sense" - " together.\n"); - return 1; - } - if ( filename == NULL ) { ERROR("Please specify filename using the -i option\n"); return 1; } + if ( output == NULL ) { + output = strdup("processed.hkl"); + } + if ( pdb == NULL ) { pdb = strdup("molecule.pdb"); } @@ -685,35 +476,6 @@ int main(int argc, char *argv[]) if ( sym == NULL ) sym = strdup("1"); - /* Show useful symmetry information */ - if ( (sym != NULL) && (in_sym != NULL) ) { - - int np = num_general_equivs(in_sym) / num_general_equivs(sym); - if ( np > 1 ) { - - STATUS("Resolving point group %s into %s " - "(%i possibilities)\n", - in_sym, sym, np); - /* Get the list of twin/Bijvoet possibilities */ - twins = get_twin_possibilities(in_sym, sym); - STATUS("Twin operation indices from %s are:", in_sym); - for ( i=0; iop); - } - STATUS("\n"); - - } else { - STATUS("No resolution required to get from %s to %s\n", - in_sym, sym); - twins = NULL; - } - - } else { - STATUS("No twin resolution requested (use -a otherwise).\n"); - twins = NULL; - in_sym = strdup(sym); - } - /* Open the data stream */ if ( strcmp(filename, "-") == 0 ) { fh = stdin; @@ -726,33 +488,13 @@ int main(int argc, char *argv[]) return 1; } - /* Read the reference reflections */ - if ( reference != NULL ) { - reference_i = new_list_intensity(); - reference_items = read_reflections(reference, reference_i, - NULL, NULL, NULL); - if ( reference_items == NULL ) { - ERROR("Couldn't read '%s'\n", reference); - return 1; - } - } else { - reference_items = NULL; - reference_i = NULL; - } - - if ( outstream != NULL ) { - outfh = fopen(outstream, "w"); - if ( outfh == NULL ) { - ERROR("Couldn't open '%s'\n", outstream); - return 1; - } - } - /* Count the number of patterns in the file */ n_total_patterns = count_patterns(fh); STATUS("There are %i patterns to process\n", n_total_patterns); rewind(fh); + model = reflist_new(); + if ( histo != NULL ) { int r; r = sscanf(histo, "%i,%i,%i", &hist_h, &hist_k, &hist_l); @@ -771,13 +513,10 @@ int main(int argc, char *argv[]) } hist_i = 0; - merge_all(fh, &model, &observed, &counts, - config_maxonly, config_scale, config_sum, + merge_all(fh, model, config_maxonly, config_scale, config_sum, config_startafter, config_stopafter, - twins, in_sym, sym, n_total_patterns, - reference_items, reference_i, - hist_vals, hist_h, hist_k, hist_l, &hist_i, NULL, NULL, - outfh); + sym, n_total_patterns, + hist_vals, hist_h, hist_k, hist_l, &hist_i); rewind(fh); if ( space_for_hist && (hist_i >= space_for_hist) ) { ERROR("Histogram array was too small!\n"); @@ -790,60 +529,17 @@ int main(int argc, char *argv[]) } STATUS("Extra pass to calculate ESDs...\n"); - devs = new_list_intensity(); - esds = new_list_sigma(); rewind(fh); - merge_all(fh, &model, &observed, &counts, - config_maxonly, config_scale, 0, - config_startafter, config_stopafter, - twins, in_sym, sym, - n_total_patterns, reference_items, reference_i, - NULL, 0, 0, 0, NULL, devs, model, NULL); + merge_all(fh, model, config_maxonly, config_scale, 0, + config_startafter, config_stopafter, sym, n_total_patterns, + NULL, 0, 0, 0, NULL); - for ( i=0; ih; k = it->k, l = it->l; - - dev = lookup_intensity(devs, h, k, l); - count = lookup_count(counts, h, k, l); - - if ( count < 2 ) { - /* If we have only one measurement, the error is 100% */ - esd = lookup_sigma(model, h, k, l); - } else { - esd = sqrt(dev) / (double)count; - } - set_sigma(esds, h, k, l, esd); - - } - - if ( output != NULL ) { - - double adu_per_photon; - - if ( beam == NULL ) { - adu_per_photon = 167.0; - STATUS("No beam parameters file provided (use -b), " - "so I'm assuming 167.0 ADU per photon.\n"); - } else { - adu_per_photon = beam->adu_per_photon; - } - - write_reflections(output, observed, model, esds, - NULL, counts, cell); - - } + write_reflist(output, model, cell); + cell_free(cell); fclose(fh); free(sym); - free(model); - free(counts); + reflist_free(model); free(output); if ( cell != NULL ) cell_free(cell); diff --git a/src/reflist-utils.c b/src/reflist-utils.c new file mode 100644 index 00000000..33d373b8 --- /dev/null +++ b/src/reflist-utils.c @@ -0,0 +1,72 @@ +/* + * reflist-utils.c + * + * Utilities to complement the core reflist.c + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + + +#include + + +#include "reflist.h" +#include "cell.h" + + +void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell) +{ + Reflection *refl; + RefListIterator *iter; + + fprintf(fh, " h k l I phase sigma(I) " + " 1/d(nm^-1) counts fs/px ss/px\n"); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + signed int h, k, l; + double intensity, esd_i, s; + double fs, ss; + + get_indices(refl, &h, &k, &l); + get_detector_pos(refl, &fs, &ss); + intensity = get_intensity(refl); + esd_i = 0.0; /* FIXME! */ + s = 0.0; /* FIXME! */ + + /* h, k, l, I, sigma(I), s */ + fprintf(fh, + "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i %6.1f %6.1f\n", + h, k, l, intensity, " -", esd_i, s/1.0e9, 1, + fs, ss); + + } +} + + +int write_reflist(const char *filename, RefList *list, UnitCell *cell) +{ + FILE *fh; + + if ( filename == NULL ) { + fh = stdout; + } else { + fh = fopen(filename, "w"); + } + + if ( fh == NULL ) { + ERROR("Couldn't open output file '%s'.\n", filename); + return 1; + } + + write_reflections_to_file(fh, list, cell); + + fclose(fh); + + return 0; +} diff --git a/src/reflist-utils.h b/src/reflist-utils.h new file mode 100644 index 00000000..d5c1e655 --- /dev/null +++ b/src/reflist-utils.h @@ -0,0 +1,30 @@ +/* + * reflist-utils.h + * + * Utilities to complement the core reflist.c + * + * (c) 2006-2011 Thomas White + * + * Part of CrystFEL - crystallography with a FEL + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#ifndef REFLIST_UTILS_H +#define REFLIST_UTILS_H + + +#include "reflist.h" +#include "cell.h" + + +extern void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell); + +extern int write_reflist(const char *filename, RefList *list, + UnitCell *cell); + + +#endif /* REFLIST_UTILS_H */ diff --git a/src/reflist.c b/src/reflist.c index 3a53330c..3181b418 100644 --- a/src/reflist.c +++ b/src/reflist.c @@ -45,6 +45,13 @@ struct _refldata { /* Intensity */ double intensity; double esd_i; + + /* Redundancy */ + int redundancy; + + /* Total squared difference between all estimates of this reflection + * and the estimated mean value */ + double sum_squared_dev; }; @@ -225,6 +232,24 @@ int get_scalable(Reflection *refl) } +int get_redundancy(Reflection *refl) +{ + return refl->data.redundancy; +} + + +double get_sum_squared_dev(Reflection *refl) +{ + return refl->data.sum_squared_dev; +} + + +double get_esd_intensity(Reflection *refl) +{ + return refl->data.esd_i; +} + + /********************************** Setters ***********************************/ void set_detector_pos(Reflection *refl, double exerr, double x, double y) @@ -258,6 +283,24 @@ void set_scalable(Reflection *refl, int scalable) } +void set_redundancy(Reflection *refl, int red) +{ + refl->data.redundancy = red; +} + + +void set_sum_squared_dev(Reflection *refl, double dev) +{ + refl->data.sum_squared_dev = dev; +} + + +void set_esd_intensity(Reflection *refl, double esd) +{ + refl->data.esd_i = esd; +} + + /********************************* Insertion **********************************/ static void insert_node(Reflection *head, Reflection *new) diff --git a/src/reflist.h b/src/reflist.h index e25a16ad..f6d01ecf 100644 --- a/src/reflist.h +++ b/src/reflist.h @@ -41,6 +41,9 @@ extern double get_intensity(Reflection *refl); extern void get_partial(Reflection *refl, double *r1, double *r2, double *p, int *clamp_low, int *clamp_high); extern int get_scalable(Reflection *refl); +extern int get_redundancy(Reflection *refl); +extern double get_sum_squared_dev(Reflection *refl); +extern double get_esd_intensity(Reflection *refl); /* Set */ extern void set_detector_pos(Reflection *refl, double exerr, @@ -49,6 +52,9 @@ extern void set_partial(Reflection *refl, double r1, double r2, double p, double clamp_low, double clamp_high); extern void set_int(Reflection *refl, double intensity); extern void set_scalable(Reflection *refl, int scalable); +extern void set_redundancy(Reflection *refl, int red); +extern void set_sum_squared_dev(Reflection *refl, double dev); +extern void set_esd_intensity(Reflection *refl, double esd); /* Insertion */ extern Reflection *add_refl(RefList *list, INDICES); diff --git a/src/stream.c b/src/stream.c index 2a75216d..4b2bea9f 100644 --- a/src/stream.c +++ b/src/stream.c @@ -23,6 +23,7 @@ #include "image.h" #include "stream.h" #include "reflist.h" +#include "reflist-utils.h" #define CHUNK_START_MARKER "----- Begin chunk -----" @@ -32,7 +33,6 @@ #define REFLECTION_START_MARKER "Reflections measured after indexing" #define REFLECTION_END_MARKER "End of reflections" - static void exclusive(const char *a, const char *b) { ERROR("The stream options '%s' and '%s' are mutually exclusive.\n", @@ -103,10 +103,9 @@ int count_patterns(FILE *fh) rval = fgets(line, 1023, fh); if ( rval == NULL ) continue; - if ( (strncmp(line, "Reflections from indexing", 25) == 0) - || (strncmp(line, "New pattern", 11) == 0) ) { - n_total_patterns++; - } + chomp(line); + if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_total_patterns++; + } while ( rval != NULL ); return n_total_patterns; @@ -147,6 +146,7 @@ static UnitCell *read_orientation_matrix(FILE *fh) static int read_reflections(FILE *fh, struct image *image) { char *rval = NULL; + int first = 1; image->reflections = reflist_new(); @@ -164,17 +164,20 @@ static int read_reflections(FILE *fh, struct image *image) if ( rval == NULL ) continue; chomp(line); - if ( strcmp(line, PEAK_LIST_END_MARKER) == 0 ) return 0; + if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return 0; r = sscanf(line, "%i %i %i %f %s %f %f %i %f %f", &h, &k, &l, &intensity, phs, &sigma, &res, &cts, &fs, &ss); - if ( r != 10 ) return 1; - - refl = add_refl(image->reflections, h, k, l); - set_int(refl, intensity); - set_detector_pos(refl, fs, ss, 0.0); - /* FIXME: Set ESD */ + if ( (r != 10) && (!first) ) return 1; + + first = 0; + if ( r == 10 ) { + refl = add_refl(image->reflections, h, k, l); + set_int(refl, intensity); + set_detector_pos(refl, fs, ss, 0.0); + set_esd_intensity(refl, sigma); + } } while ( rval != NULL ); @@ -184,45 +187,10 @@ static int read_reflections(FILE *fh, struct image *image) } -static void write_reflections(struct image *image, FILE *ofh) -{ - Reflection *refl; - RefListIterator *iter; - - fprintf(ofh, REFLECTION_START_MARKER"\n"); - /* FIXME: Unify this with write_reflections() over in reflections.c */ - fprintf(ofh, " h k l I phase sigma(I) " - " 1/d(nm^-1) counts fs/px ss/px\n"); - - for ( refl = first_refl(image->reflections, &iter); - refl != NULL; - refl = next_refl(refl, iter) ) { - - signed int h, k, l; - double intensity, esd_i, s; - double fs, ss; - - get_indices(refl, &h, &k, &l); - get_detector_pos(refl, &fs, &ss); - intensity = get_intensity(refl); - esd_i = 0.0; /* FIXME! */ - s = 0.0; /* FIXME! */ - - /* h, k, l, I, sigma(I), s */ - fprintf(ofh, - "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i %6.1f %6.1f\n", - h, k, l, intensity, " -", esd_i, s/1.0e9, 1, - fs, ss); - - } - - fprintf(ofh, REFLECTION_END_MARKER"\n"); -} - - static int read_peaks(FILE *fh, struct image *image) { char *rval = NULL; + int first = 1; image->features = image_feature_list_new(); @@ -239,9 +207,17 @@ static int read_peaks(FILE *fh, struct image *image) if ( strcmp(line, PEAK_LIST_END_MARKER) == 0 ) return 0; r = sscanf(line, "%f %f %f %f", &x, &y, &d, &intensity); - if ( r != 4 ) return 1; + if ( (r != 4) && (!first) ) { + ERROR("Failed to parse peak list line.\n"); + ERROR("The failed line was: '%s'\n", line); + return 1; + } - image_add_feature(image->features, x, y, image, 1.0, NULL); + first = 0; + if ( r == 4 ) { + image_add_feature(image->features, x, y, + image, 1.0, NULL); + } } while ( rval != NULL ); @@ -255,7 +231,7 @@ static void write_peaks(struct image *image, FILE *ofh) int i; fprintf(ofh, PEAK_LIST_START_MARKER"\n"); - fprintf(ofh, " fs/px ss/px (1/d)/nm^-1 Intensity\n"); + fprintf(ofh, " fs/px ss/px (1/d)/nm^-1 Intensity\n"); for ( i=0; ifeatures); i++ ) { @@ -269,7 +245,7 @@ static void write_peaks(struct image *image, FILE *ofh) r = get_q(image, f->fs, f->ss, NULL, 1.0/image->lambda); q = modulus(r.u, r.v, r.w); - fprintf(ofh, "%8.3f %8.3f %8.3f %12.3f\n", + fprintf(ofh, "%6.1f %6.1f %10.2f %10.2f\n", f->fs, f->ss, q/1.0e9, f->intensity); } @@ -330,8 +306,12 @@ void write_chunk(FILE *ofh, struct image *i, int f) } if ( (f & STREAM_PIXELS) || (f & STREAM_INTEGRATED) ) { + fprintf(ofh, "\n"); - write_reflections(i, ofh); + fprintf(ofh, REFLECTION_START_MARKER"\n"); + write_reflections_to_file(ofh, i->reflections, i->indexed_cell); + fprintf(ofh, REFLECTION_END_MARKER"\n"); + } fprintf(ofh, CHUNK_END_MARKER"\n\n"); @@ -374,10 +354,11 @@ int read_chunk(FILE *fh, struct image *image) if ( find_start_of_chunk(fh) ) return 1; image->i0_available = 0; - if ( image->features != NULL ) { - image_feature_list_free(image->features); - image->features = NULL; - } + image->i0 = 1.0; + image->lambda = -1.0; + image->features = NULL; + image->reflections = NULL; + image->indexed_cell = NULL; do { @@ -426,16 +407,22 @@ int read_chunk(FILE *fh, struct image *image) } if ( strcmp(line, PEAK_LIST_START_MARKER) == 0 ) { - if ( read_peaks(fh, image) ) return 1; + if ( read_peaks(fh, image) ) { + ERROR("Failed while reading peaks\n"); + return 1; + } } if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) { - if ( read_reflections(fh, image) ) return 1; + if ( read_reflections(fh, image) ) { + ERROR("Failed while reading reflections\n"); + return 1; + } } } while ( strcmp(line, CHUNK_END_MARKER) != 0 ); - if ( have_filename && have_cell && have_ev ) return 0; + if ( have_filename && have_ev ) return 0; ERROR("Incomplete chunk found in input file.\n"); return 1; @@ -496,3 +483,24 @@ int find_chunk(FILE *fh, UnitCell **cell, char **filename, double *ev) return 1; } + + +int skip_some_files(FILE *fh, int n) +{ + char *rval = NULL; + int n_patterns = 0; + + do { + + char line[1024]; + + if ( n_patterns == n ) return 0; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + if ( strcmp(line, CHUNK_END_MARKER) == 0 ) n_patterns++; + + } while ( rval != NULL ); + + return 1; +} diff --git a/src/stream.h b/src/stream.h index 3869e515..64bd9529 100644 --- a/src/stream.h +++ b/src/stream.h @@ -38,4 +38,8 @@ extern void write_chunk(FILE *ofh, struct image *image, int flags); extern int parse_stream_flags(const char *a); +extern int read_chunk(FILE *fh, struct image *image); + +extern int skip_some_files(FILE *fh, int n); + #endif /* STREAM_H */ -- cgit v1.2.3