diff options
| -rw-r--r-- | Makefile.am | 32 | ||||
| -rw-r--r-- | Makefile.in | 76 | ||||
| -rw-r--r-- | src/check_hkl.c | 92 | ||||
| -rw-r--r-- | src/compare_hkl.c | 220 | ||||
| -rw-r--r-- | src/get_hkl.c | 283 | ||||
| -rw-r--r-- | src/indexamajig.c | 2 | ||||
| -rw-r--r-- | src/partialator.c | 4 | ||||
| -rw-r--r-- | src/pattern_sim.c | 30 | ||||
| -rw-r--r-- | src/povray.c | 47 | ||||
| -rw-r--r-- | src/povray.h | 8 | ||||
| -rw-r--r-- | src/process_hkl.c | 9 | ||||
| -rw-r--r-- | src/reflections.c | 230 | ||||
| -rw-r--r-- | src/reflections.h | 38 | ||||
| -rw-r--r-- | src/reflist-utils.c | 247 | ||||
| -rw-r--r-- | src/reflist-utils.h | 18 | ||||
| -rw-r--r-- | src/reflist.c | 27 | ||||
| -rw-r--r-- | src/reflist.h | 6 | ||||
| -rw-r--r-- | src/render_hkl.c | 53 | ||||
| -rw-r--r-- | src/statistics.c | 385 | ||||
| -rw-r--r-- | src/statistics.h | 48 | ||||
| -rw-r--r-- | src/stream.c | 49 |
21 files changed, 979 insertions, 925 deletions
diff --git a/Makefile.am b/Makefile.am index 50631ced..2173bfe4 100644 --- a/Makefile.am +++ b/Makefile.am @@ -26,7 +26,7 @@ LDADD = $(top_builddir)/lib/libgnu.a @IGNORE_UNUSED_LIBRARIES_CFLAGS@ src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c src/utils.c \ src/image.c src/cell.c src/hdf5-file.c \ src/detector.c src/peaks.c \ - src/reflections.c src/beam-parameters.c \ + src/reflist-utils.c src/beam-parameters.c \ src/symmetry.c src/thread-pool.c src/reflist.c if HAVE_OPENCL src_pattern_sim_SOURCES += src/diffraction-gpu.c src/cl-utils.c @@ -51,10 +51,10 @@ 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/dirax.c src/mosflm.c \ - src/reflections.c src/symmetry.c \ + src/reflist-utils.c src/symmetry.c \ src/geometry.c src/thread-pool.c \ - src/beam-parameters.c src/reflist.c src/stream.c \ - src/reflist-utils.c + src/beam-parameters.c src/reflist.c src/stream.c + if HAVE_OPENCL src_indexamajig_SOURCES += src/diffraction-gpu.c src/cl-utils.c endif @@ -67,25 +67,25 @@ src_hdfsee_SOURCES = src/hdfsee.c src/dw-hdfsee.c src/render.c \ endif src_get_hkl_SOURCES = src/get_hkl.c src/cell.c src/utils.c \ - src/reflections.c src/symmetry.c src/beam-parameters.c \ - src/thread-pool.c + src/reflist-utils.c src/symmetry.c src/beam-parameters.c \ + src/thread-pool.c src/reflist.c src_compare_hkl_SOURCES = src/compare_hkl.c src/cell.c src/utils.c \ - src/reflections.c src/statistics.c src/symmetry.c \ - src/thread-pool.c + src/reflist-utils.c src/statistics.c src/symmetry.c \ + src/thread-pool.c src/reflist.c src_check_hkl_SOURCES = src/check_hkl.c src/cell.c src/utils.c \ - src/reflections.c src/statistics.c src/symmetry.c \ - src/thread-pool.c + src/reflist-utils.c src/statistics.c src/symmetry.c \ + src/thread-pool.c src/reflist.c src_powder_plot_SOURCES = src/powder_plot.c src/cell.c src/utils.c src/image.c \ src/hdf5-file.c src/detector.c src/thread-pool.c \ src/beam-parameters.c -src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflections.c \ - src/utils.c src/povray.c src/symmetry.c src/render.c \ +src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflist-utils.c \ + src/utils.c src/povray.c src/symmetry.c src/render.c \ src/hdf5-file.c src/image.c src/filters.c \ - src/thread-pool.c + src/thread-pool.c src/reflist.c src_sum_stack_SOURCES = src/sum_stack.c src/utils.c src/hdf5-file.c \ src/image.c src/filters.c src/peaks.c src/detector.c \ @@ -97,10 +97,10 @@ src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \ src/utils.c src/detector.c src/peaks.c src/image.c \ - src/geometry.c src/reflections.c src/stream.c \ + src/geometry.c src/reflist-utils.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/reflist-utils.c + src/hrs-scaling.c src/reflist.c if BUILD_CUBEIT src_cubeit_SOURCES = src/cubeit.c src/cell.c src/hdf5-file.c src/utils.c \ @@ -115,7 +115,7 @@ tests_list_check_SOURCES = tests/list_check.c src/reflist.c src/thread-pool.c \ INCLUDES = "-I$(top_srcdir)/data" EXTRA_DIST += src/cell.h src/hdf5-file.h src/image.h src/utils.h \ - src/diffraction.h src/detector.h src/reflections.h \ + src/diffraction.h src/detector.h src/reflist-utils.h \ src/list_tmp.h src/statistics.h src/dw-hdfsee.h \ src/render.h src/hdfsee.h src/dirax.h src/mosflm.h src/peaks.h \ src/index.h src/filters.h src/diffraction-gpu.h src/cl-utils.h \ diff --git a/Makefile.in b/Makefile.in index f8d4bf23..226389f9 100644 --- a/Makefile.in +++ b/Makefile.in @@ -90,16 +90,17 @@ src_calibrate_detector_OBJECTS = $(am_src_calibrate_detector_OBJECTS) src_calibrate_detector_LDADD = $(LDADD) src_calibrate_detector_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_src_check_hkl_OBJECTS = src/check_hkl.$(OBJEXT) src/cell.$(OBJEXT) \ - src/utils.$(OBJEXT) src/reflections.$(OBJEXT) \ + src/utils.$(OBJEXT) src/reflist-utils.$(OBJEXT) \ src/statistics.$(OBJEXT) src/symmetry.$(OBJEXT) \ - src/thread-pool.$(OBJEXT) + src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT) src_check_hkl_OBJECTS = $(am_src_check_hkl_OBJECTS) src_check_hkl_LDADD = $(LDADD) src_check_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_src_compare_hkl_OBJECTS = src/compare_hkl.$(OBJEXT) \ src/cell.$(OBJEXT) src/utils.$(OBJEXT) \ - src/reflections.$(OBJEXT) src/statistics.$(OBJEXT) \ - src/symmetry.$(OBJEXT) src/thread-pool.$(OBJEXT) + src/reflist-utils.$(OBJEXT) src/statistics.$(OBJEXT) \ + src/symmetry.$(OBJEXT) src/thread-pool.$(OBJEXT) \ + src/reflist.$(OBJEXT) src_compare_hkl_OBJECTS = $(am_src_compare_hkl_OBJECTS) src_compare_hkl_LDADD = $(LDADD) src_compare_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -120,9 +121,9 @@ src_cubeit_OBJECTS = $(am_src_cubeit_OBJECTS) src_cubeit_LDADD = $(LDADD) src_cubeit_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_src_get_hkl_OBJECTS = src/get_hkl.$(OBJEXT) src/cell.$(OBJEXT) \ - src/utils.$(OBJEXT) src/reflections.$(OBJEXT) \ + src/utils.$(OBJEXT) src/reflist-utils.$(OBJEXT) \ src/symmetry.$(OBJEXT) src/beam-parameters.$(OBJEXT) \ - src/thread-pool.$(OBJEXT) + src/thread-pool.$(OBJEXT) src/reflist.$(OBJEXT) src_get_hkl_OBJECTS = $(am_src_get_hkl_OBJECTS) src_get_hkl_LDADD = $(LDADD) src_get_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -142,10 +143,9 @@ src_hdfsee_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am__src_indexamajig_SOURCES_DIST = 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/dirax.c \ - src/mosflm.c src/reflections.c src/symmetry.c src/geometry.c \ + src/mosflm.c src/reflist-utils.c src/symmetry.c src/geometry.c \ src/thread-pool.c src/beam-parameters.c src/reflist.c \ - src/stream.c src/reflist-utils.c src/diffraction-gpu.c \ - src/cl-utils.c + src/stream.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) \ @@ -153,35 +153,34 @@ am_src_indexamajig_OBJECTS = src/indexamajig.$(OBJEXT) \ src/image.$(OBJEXT) src/peaks.$(OBJEXT) src/index.$(OBJEXT) \ src/filters.$(OBJEXT) src/diffraction.$(OBJEXT) \ src/detector.$(OBJEXT) src/dirax.$(OBJEXT) \ - src/mosflm.$(OBJEXT) src/reflections.$(OBJEXT) \ + src/mosflm.$(OBJEXT) src/reflist-utils.$(OBJEXT) \ src/symmetry.$(OBJEXT) src/geometry.$(OBJEXT) \ src/thread-pool.$(OBJEXT) src/beam-parameters.$(OBJEXT) \ - src/reflist.$(OBJEXT) src/stream.$(OBJEXT) \ - src/reflist-utils.$(OBJEXT) $(am__objects_1) + src/reflist.$(OBJEXT) src/stream.$(OBJEXT) $(am__objects_1) src_indexamajig_OBJECTS = $(am_src_indexamajig_OBJECTS) src_indexamajig_LDADD = $(LDADD) src_indexamajig_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_src_partialator_OBJECTS = src/partialator.$(OBJEXT) \ src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) src/utils.$(OBJEXT) \ src/detector.$(OBJEXT) src/peaks.$(OBJEXT) src/image.$(OBJEXT) \ - src/geometry.$(OBJEXT) src/reflections.$(OBJEXT) \ + src/geometry.$(OBJEXT) src/reflist-utils.$(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-utils.$(OBJEXT) + src/reflist.$(OBJEXT) src_partialator_OBJECTS = $(am_src_partialator_OBJECTS) src_partialator_LDADD = $(LDADD) src_partialator_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am__src_pattern_sim_SOURCES_DIST = src/pattern_sim.c src/diffraction.c \ src/utils.c src/image.c src/cell.c src/hdf5-file.c \ - src/detector.c src/peaks.c src/reflections.c \ + src/detector.c src/peaks.c src/reflist-utils.c \ src/beam-parameters.c src/symmetry.c src/thread-pool.c \ src/reflist.c src/diffraction-gpu.c src/cl-utils.c am_src_pattern_sim_OBJECTS = src/pattern_sim.$(OBJEXT) \ src/diffraction.$(OBJEXT) src/utils.$(OBJEXT) \ src/image.$(OBJEXT) src/cell.$(OBJEXT) src/hdf5-file.$(OBJEXT) \ src/detector.$(OBJEXT) src/peaks.$(OBJEXT) \ - src/reflections.$(OBJEXT) src/beam-parameters.$(OBJEXT) \ + src/reflist-utils.$(OBJEXT) src/beam-parameters.$(OBJEXT) \ src/symmetry.$(OBJEXT) src/thread-pool.$(OBJEXT) \ src/reflist.$(OBJEXT) $(am__objects_1) src_pattern_sim_OBJECTS = $(am_src_pattern_sim_OBJECTS) @@ -205,11 +204,12 @@ src_process_hkl_OBJECTS = $(am_src_process_hkl_OBJECTS) src_process_hkl_LDADD = $(LDADD) src_process_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a am_src_render_hkl_OBJECTS = src/render_hkl.$(OBJEXT) \ - src/cell.$(OBJEXT) src/reflections.$(OBJEXT) \ + src/cell.$(OBJEXT) src/reflist-utils.$(OBJEXT) \ src/utils.$(OBJEXT) src/povray.$(OBJEXT) \ src/symmetry.$(OBJEXT) src/render.$(OBJEXT) \ src/hdf5-file.$(OBJEXT) src/image.$(OBJEXT) \ - src/filters.$(OBJEXT) src/thread-pool.$(OBJEXT) + src/filters.$(OBJEXT) src/thread-pool.$(OBJEXT) \ + src/reflist.$(OBJEXT) src_render_hkl_OBJECTS = $(am_src_render_hkl_OBJECTS) src_render_hkl_LDADD = $(LDADD) src_render_hkl_DEPENDENCIES = $(top_builddir)/lib/libgnu.a @@ -558,7 +558,7 @@ top_builddir = @top_builddir@ top_srcdir = @top_srcdir@ EXTRA_DIST = configure m4/gnulib-cache.m4 src/cell.h src/hdf5-file.h \ src/image.h src/utils.h src/diffraction.h src/detector.h \ - src/reflections.h src/list_tmp.h src/statistics.h \ + src/reflist-utils.h src/list_tmp.h src/statistics.h \ src/dw-hdfsee.h src/render.h src/hdfsee.h src/dirax.h \ src/mosflm.h src/peaks.h src/index.h src/filters.h \ src/diffraction-gpu.h src/cl-utils.h src/symmetry.h \ @@ -588,7 +588,7 @@ AM_CPPFLAGS = -DDATADIR=\""$(datadir)"\" -I$(top_builddir)/lib -I$(top_srcdir)/l LDADD = $(top_builddir)/lib/libgnu.a @IGNORE_UNUSED_LIBRARIES_CFLAGS@ src_pattern_sim_SOURCES = src/pattern_sim.c src/diffraction.c \ src/utils.c src/image.c src/cell.c src/hdf5-file.c \ - src/detector.c src/peaks.c src/reflections.c \ + src/detector.c src/peaks.c src/reflist-utils.c \ src/beam-parameters.c src/symmetry.c src/thread-pool.c \ src/reflist.c $(am__append_3) @HAVE_OPENCL_TRUE@tests_gpu_sim_check_SOURCES = tests/gpu_sim_check.c src/utils.c \ @@ -606,33 +606,33 @@ src_process_hkl_SOURCES = src/process_hkl.c src/statistics.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/dirax.c \ - src/mosflm.c src/reflections.c src/symmetry.c src/geometry.c \ + src/mosflm.c src/reflist-utils.c src/symmetry.c src/geometry.c \ src/thread-pool.c src/beam-parameters.c src/reflist.c \ - src/stream.c src/reflist-utils.c $(am__append_6) + src/stream.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 src_get_hkl_SOURCES = src/get_hkl.c src/cell.c src/utils.c \ - src/reflections.c src/symmetry.c src/beam-parameters.c \ - src/thread-pool.c + src/reflist-utils.c src/symmetry.c src/beam-parameters.c \ + src/thread-pool.c src/reflist.c src_compare_hkl_SOURCES = src/compare_hkl.c src/cell.c src/utils.c \ - src/reflections.c src/statistics.c src/symmetry.c \ - src/thread-pool.c + src/reflist-utils.c src/statistics.c src/symmetry.c \ + src/thread-pool.c src/reflist.c src_check_hkl_SOURCES = src/check_hkl.c src/cell.c src/utils.c \ - src/reflections.c src/statistics.c src/symmetry.c \ - src/thread-pool.c + src/reflist-utils.c src/statistics.c src/symmetry.c \ + src/thread-pool.c src/reflist.c src_powder_plot_SOURCES = src/powder_plot.c src/cell.c src/utils.c src/image.c \ src/hdf5-file.c src/detector.c src/thread-pool.c \ src/beam-parameters.c -src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflections.c \ - src/utils.c src/povray.c src/symmetry.c src/render.c \ +src_render_hkl_SOURCES = src/render_hkl.c src/cell.c src/reflist-utils.c \ + src/utils.c src/povray.c src/symmetry.c src/render.c \ src/hdf5-file.c src/image.c src/filters.c \ - src/thread-pool.c + src/thread-pool.c src/reflist.c src_sum_stack_SOURCES = src/sum_stack.c src/utils.c src/hdf5-file.c \ src/image.c src/filters.c src/peaks.c src/detector.c \ @@ -644,10 +644,10 @@ src_calibrate_detector_SOURCES = src/calibrate_detector.c src/utils.c \ src_partialator_SOURCES = src/partialator.c src/cell.c src/hdf5-file.c \ src/utils.c src/detector.c src/peaks.c src/image.c \ - src/geometry.c src/reflections.c src/stream.c \ + src/geometry.c src/reflist-utils.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/reflist-utils.c + src/hrs-scaling.c src/reflist.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 \ @@ -803,12 +803,14 @@ src/calibrate_detector$(EXEEXT): $(src_calibrate_detector_OBJECTS) $(src_calibra src/check_hkl.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/cell.$(OBJEXT): src/$(am__dirstamp) src/$(DEPDIR)/$(am__dirstamp) -src/reflections.$(OBJEXT): src/$(am__dirstamp) \ +src/reflist-utils.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/statistics.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) src/symmetry.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) +src/reflist.$(OBJEXT): src/$(am__dirstamp) \ + src/$(DEPDIR)/$(am__dirstamp) src/check_hkl$(EXEEXT): $(src_check_hkl_OBJECTS) $(src_check_hkl_DEPENDENCIES) src/$(am__dirstamp) @rm -f src/check_hkl$(EXEEXT) $(AM_V_CCLD)$(LINK) $(src_check_hkl_OBJECTS) $(src_check_hkl_LDADD) $(LIBS) @@ -825,10 +827,6 @@ src/filters.$(OBJEXT): src/$(am__dirstamp) \ src/$(DEPDIR)/$(am__dirstamp) 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) @@ -948,7 +946,6 @@ mostlyclean-compile: -rm -f src/povray.$(OBJEXT) -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) @@ -994,7 +991,6 @@ distclean-compile: @AMDEP_TRUE@@am__include@ @am__quote@src/$(DEPDIR)/povray.Po@am__quote@ @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@ diff --git a/src/check_hkl.c b/src/check_hkl.c index fb6d02d5..903188a1 100644 --- a/src/check_hkl.c +++ b/src/check_hkl.c @@ -22,9 +22,10 @@ #include <getopt.h> #include "utils.h" -#include "reflections.h" #include "statistics.h" #include "symmetry.h" +#include "reflist.h" +#include "reflist-utils.h" /* Number of bins for plot of resolution shells */ @@ -46,9 +47,8 @@ static void show_help(const char *s) } -static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell, - const char *sym, unsigned int *counts, - const double *sigma, double rmin_fix, double rmax_fix) +static void plot_shells(RefList *list, UnitCell *cell, const char *sym, + double rmin_fix, double rmax_fix) { double num[NBINS]; int cts[NBINS]; @@ -70,6 +70,8 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell, int nmeas = 0; int nmeastot = 0; int nout = 0; + Reflection *refl; + RefListIterator *iter; if ( cell == NULL ) { ERROR("Need the unit cell to plot resolution shells.\n"); @@ -93,16 +95,17 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell, mean[i] = 0; } - rmin = +INFINITY; - rmax = 0.0; - for ( i=0; i<num_items(items); i++ ) { + /* Iterate over all common reflections and calculate min and max + * resolution */ + rmin = +INFINITY; rmax = 0.0; + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; signed int h, k, l; double d; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl, &h, &k, &l); d = resolution(cell, h, k, l) * 2.0; if ( d > rmax ) rmax = d; @@ -179,16 +182,16 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell, free(counted); /* Calculate means */ - for ( i=0; i<num_items(items); i++ ) { + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; signed int h, k, l; double d; int bin; int j; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl, &h, &k, &l); d = resolution(cell, h, k, l) * 2.0; @@ -205,7 +208,7 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell, } measured[bin]++; - mean[bin] += lookup_intensity(ref, h, k, l); + mean[bin] += get_intensity(refl); } @@ -214,21 +217,21 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell, } /* Characterise the data set */ - for ( i=0; i<num_items(items); i++ ) { + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; signed int h, k, l; double d; int bin; int j; double val, esd; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl, &h, &k, &l); d = resolution(cell, h, k, l) * 2.0; - val = lookup_intensity(ref, h, k, l); - esd = lookup_sigma(sigma, h, k, l); + val = get_intensity(refl); + esd = get_esd_intensity(refl); bin = -1; for ( j=0; j<NBINS; j++ ) { @@ -245,11 +248,11 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell, if ( !isfinite(val/esd) ) continue; /* measured[bin] was done earlier */ - measurements[bin] += lookup_count(counts, h, k, l); + measurements[bin] += get_redundancy(refl); snr[bin] += val / esd; snr_total += val / esd; nmeas++; - nmeastot += lookup_count(counts, h, k, l); + nmeastot += get_redundancy(refl); var[bin] += pow(val-mean[bin], 2.0); @@ -289,17 +292,15 @@ static void plot_shells(const double *ref, ReflItemList *items, UnitCell *cell, int main(int argc, char *argv[]) { int c; - double *ref; UnitCell *cell; char *file = NULL; char *sym = NULL; - int i; - ReflItemList *items; - ReflItemList *good_items; + RefList *raw_list; + RefList *list; + Reflection *refl; + RefListIterator *iter; char *pdb = NULL; - double *esd; int rej = 0; - unsigned int *cts; float rmin_fix = -1.0; float rmax_fix = -1.0; @@ -369,37 +370,35 @@ int main(int argc, char *argv[]) cell = load_cell_from_pdb(pdb); free(pdb); - ref = new_list_intensity(); - esd = new_list_sigma(); - cts = new_list_count(); - items = read_reflections(file, ref, NULL, cts, esd); - if ( items == NULL ) { - ERROR("Couldn't open file '%s'\n", file); + raw_list = read_reflections(file); + if ( raw_list == NULL ) { + ERROR("Couldn't read file '%s'\n", file); return 1; } /* Check that the intensities have the correct symmetry */ - if ( check_symmetry(items, sym) ) { + if ( check_list_symmetry(raw_list, sym) ) { ERROR("The input reflection list does not appear to" " have symmetry %s\n", sym); return 1; } - /* Reject reflections */ - good_items = new_items(); - for ( i=0; i<num_items(items); i++ ) { + /* Reject some reflections */ + list = reflist_new(); + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; signed int h, k, l; double val, sig; int ig = 0; double d; + Reflection *new; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl, &h, &k, &l); - val = lookup_intensity(ref, h, k, l); - sig = lookup_sigma(esd, h, k, l); + val = get_intensity(refl); + sig = get_esd_intensity(refl); if ( val < 3.0 * sig ) { rej++; @@ -412,11 +411,12 @@ int main(int argc, char *argv[]) //if ( ig ) continue; - add_item(good_items, h, k, l); + new = add_refl(list, h, k, l); + copy_data(new, refl); } - plot_shells(ref, items, cell, sym, cts, esd, rmin_fix, rmax_fix); + plot_shells(list, cell, sym, rmin_fix, rmax_fix); return 0; } diff --git a/src/compare_hkl.c b/src/compare_hkl.c index abe8bc6a..16a744ea 100644 --- a/src/compare_hkl.c +++ b/src/compare_hkl.c @@ -20,11 +20,12 @@ #include <string.h> #include <unistd.h> #include <getopt.h> +#include <assert.h> #include "utils.h" -#include "reflections.h" #include "statistics.h" #include "symmetry.h" +#include "reflist-utils.h" /* Number of bins for plot of resolution shells */ @@ -38,7 +39,7 @@ static void show_help(const char *s) "Compare intensity lists.\n" "\n" " -h, --help Display this help message.\n" -" -o, --output=<filename> Specify output filename for correction factor.\n" +" -o, --ratio=<filename> Specify output filename for ratios.\n" " -y, --symmetry=<sym> The symmetry of both the input files.\n" " -p, --pdb=<filename> PDB file to use (default: molecule.pdb).\n" " --shells Plot the figures of merit by resolution.\n" @@ -48,9 +49,9 @@ static void show_help(const char *s) } -static void plot_shells(const double *ref1, const double *ref2, - ReflItemList *items, double scale, UnitCell *cell, - const char *sym, double rmin_fix, double rmax_fix) +static void plot_shells(RefList *list1, RefList *list2, double scale, + UnitCell *cell, const char *sym, + double rmin_fix, double rmax_fix) { double num[NBINS]; int cts[NBINS]; @@ -69,6 +70,12 @@ static void plot_shells(const double *ref1, const double *ref2, int ctot; FILE *fh; int nout = 0; + Reflection *refl1; + RefListIterator *iter; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + int hmax, kmax, lmax; if ( cell == NULL ) { ERROR("Need the unit cell to plot resolution shells.\n"); @@ -90,16 +97,20 @@ static void plot_shells(const double *ref1, const double *ref2, snr[i] = 0; } - rmin = +INFINITY; - rmax = 0.0; - for ( i=0; i<num_items(items); i++ ) { + /* Iterate over all common reflections and calculate min and max + * resolution */ + rmin = +INFINITY; rmax = 0.0; + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - struct refl_item *it; signed int h, k, l; double d; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ d = resolution(cell, h, k, l) * 2.0; if ( d > rmax ) rmax = d; @@ -117,6 +128,7 @@ static void plot_shells(const double *ref1, const double *ref2, if ( rmin_fix > 0.0 ) rmin = rmin_fix; if ( rmax_fix > 0.0 ) rmax = rmax_fix; + /* Calculate the resolution bins */ total_vol = pow(rmax, 3.0) - pow(rmin, 3.0); vol_per_shell = total_vol / NBINS; rmins[0] = rmin; @@ -145,9 +157,16 @@ static void plot_shells(const double *ref1, const double *ref2, /* Count the number of reflections possible in each shell */ counted = new_list_count(); - for ( h=-50; h<=+50; h++ ) { - for ( k=-50; k<=+50; k++ ) { - for ( l=-50; l<=+50; l++ ) { + cell_get_reciprocal(cell, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + hmax = rmax / modulus(asx, asy, asz); + kmax = rmax / modulus(bsx, bsy, bsz); + lmax = rmax / modulus(csx, csy, csz); + + for ( h=-hmax; h<hmax; h++ ) { + for ( k=-kmax; k<kmax; k++ ) { + for ( l=-lmax; l<lmax; l++ ) { double d; signed int hs, ks, ls; @@ -178,17 +197,21 @@ static void plot_shells(const double *ref1, const double *ref2, den = 0.0; ctot = 0; nout = 0; - for ( i=0; i<num_items(items); i++ ) { - struct refl_item *it; + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { + signed int h, k, l; double d; int bin; double i1, i2; int j; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ d = resolution(cell, h, k, l) * 2.0; @@ -206,8 +229,8 @@ static void plot_shells(const double *ref1, const double *ref2, continue; } - i1 = lookup_intensity(ref1, h, k, l); - i2 = lookup_intensity(ref2, h, k, l); + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); i2 *= scale; num[bin] += fabs(i1 - i2); @@ -238,33 +261,32 @@ static void plot_shells(const double *ref1, const double *ref2, int main(int argc, char *argv[]) { int c; - double *ref1; - double *ref2; - double *ref2_transformed; - double *out; UnitCell *cell; - char *outfile = NULL; + char *ratiofile = NULL; char *afile = NULL; char *bfile = NULL; char *sym = NULL; double scale, scale_r2, scale_rdig, R1, R2, R1i, Rdiff, pearson; double scale_rintint, scale_r1i, scale_r1, scale_r1fi; - int i, ncom; - ReflItemList *i1, *i2, *icommon; + int ncom; + RefList *list1; + RefList *list2; + RefList *list2_transformed; + RefList *ratio; + RefList *deleteme; int config_shells = 0; char *pdb = NULL; - double *esd1; - double *esd2; int rej1 = 0; int rej2 = 0; - unsigned int *cts1; float rmin_fix = -1.0; float rmax_fix = -1.0; + Reflection *refl1; + RefListIterator *iter; /* Long options */ const struct option longopts[] = { {"help", 0, NULL, 'h'}, - {"output", 1, NULL, 'o'}, + {"ratio" , 1, NULL, 'o'}, {"symmetry", 1, NULL, 'y'}, {"shells", 0, &config_shells, 1}, {"pdb", 1, NULL, 'p'}, @@ -282,7 +304,7 @@ int main(int argc, char *argv[]) return 0; case 'o' : - outfile = strdup(optarg); + ratiofile = strdup(optarg); break; case 'y' : @@ -335,63 +357,62 @@ int main(int argc, char *argv[]) cell = load_cell_from_pdb(pdb); free(pdb); - ref1 = new_list_intensity(); - esd1 = new_list_sigma(); - cts1 = new_list_count(); - i1 = read_reflections(afile, ref1, NULL, cts1, esd1); - if ( i1 == NULL ) { - ERROR("Couldn't open file '%s'\n", afile); + list1 = read_reflections(afile); + if ( list1 == NULL ) { + ERROR("Couldn't read file '%s'\n", afile); return 1; } - ref2 = new_list_intensity(); - esd2 = new_list_sigma(); - i2 = read_reflections(bfile, ref2, NULL, NULL, esd2); - if ( i2 == NULL ) { - ERROR("Couldn't open file '%s'\n", bfile); + + list2 = read_reflections(bfile); + if ( list2 == NULL ) { + ERROR("Couldn't read file '%s'\n", bfile); return 1; } /* Check that the intensities have the correct symmetry */ - if ( check_symmetry(i1, sym) ) { + if ( check_list_symmetry(list1, sym) ) { ERROR("The first input reflection list does not appear to" " have symmetry %s\n", sym); return 1; } - if ( check_symmetry(i2, sym) ) { + if ( check_list_symmetry(list2, sym) ) { ERROR("The second input reflection list does not appear to" " have symmetry %s\n", sym); return 1; } - /* List for output scale factor map */ - out = new_list_intensity(); - /* Find common reflections (taking symmetry into account) */ - icommon = new_items(); - ref2_transformed = new_list_intensity(); - for ( i=0; i<num_items(i1); i++ ) { + list2_transformed = reflist_new(); + ratio = reflist_new(); + deleteme = reflist_new(); + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - struct refl_item *it; signed int h, k, l; signed int he, ke, le; double val1, val2; double sig1, sig2; int ig = 0; double d; + Reflection *refl2; + Reflection *tr; - it = get_item(i1, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); - if ( !find_unique_equiv(i2, h, k, l, sym, &he, &ke, &le) ) { - //STATUS("%i %i %i not matched (%f nm).\n", h, k, l, - // 1.0/(2.0*resolution(cell, h, k, l)/1e9)); + if ( !find_equiv_in_list(list2, h, k, l, sym, &he, &ke, &le) ) { + /* No common reflection */ + add_refl(deleteme, h, k, l); continue; } - val1 = lookup_intensity(ref1, h, k, l); - val2 = lookup_intensity(ref2, he, ke, le); - sig1 = lookup_sigma(esd1, h, k, l); - sig2 = lookup_sigma(esd2, he, ke, le); + refl2 = find_refl(list2, he, ke, le); + assert(refl2 != NULL); + + val1 = get_intensity(refl1); + val2 = get_intensity(refl2); + sig1 = get_esd_intensity(refl1); + sig2 = get_esd_intensity(refl2); if ( val1 < 3.0 * sig1 ) { rej1++; @@ -408,67 +429,90 @@ int main(int argc, char *argv[]) //if ( ig ) continue; - set_intensity(ref2_transformed, h, k, l, val2); - set_intensity(out, h, k, l, val1/val2); - add_item(icommon, h, k, l); + /* Add the old data from 'refl2' to a new list with the same + * indices as its equivalent in 'list1' */ + tr = add_refl(list2_transformed, h, k, l); + copy_data(tr, refl2); + + /* Add divided version to 'output' list */ + tr = add_refl(ratio, h, k, l); + set_int(tr, val1/val2); } - ncom = num_items(icommon); + STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej1, afile); STATUS("%i reflections in '%s' had I < 3.0*sigma(I)\n", rej2, bfile); + ncom = num_reflections(list2_transformed); STATUS("%i,%i reflections: %i in common\n", - num_items(i1), num_items(i2), ncom); + num_reflections(list1), num_reflections(list2), ncom); + + /* Trim reflections from 'list1' which had no equivalents in 'list2' */ + for ( refl1 = first_refl(deleteme, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - R1 = stat_r1_ignore(ref1, ref2_transformed, icommon, &scale_r1fi); - STATUS("R1(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n", + signed int h, k, l; + Reflection *del; + + get_indices(refl1, &h, &k, &l); + del = find_refl(list1, h, k, l); + assert(del != NULL); + + delete_refl(del); + + } + reflist_free(deleteme); + reflist_free(list2); + + R1 = stat_r1_ignore(list1, list2_transformed, &scale_r1fi); + STATUS("R1(F) = %5.4f %% (scale=%5.2e)" + " (ignoring negative intensities)\n", R1*100.0, scale_r1fi); - R1 = stat_r1_zero(ref1, ref2_transformed, icommon, &scale_r1); - STATUS("R1(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n", + R1 = stat_r1_zero(list1, list2_transformed, &scale_r1); + STATUS("R1(F) = %5.4f %% (scale=%5.2e)" + " (zeroing negative intensities)\n", R1*100.0, scale_r1); - R2 = stat_r2(ref1, ref2_transformed, icommon, &scale_r2); + R2 = stat_r2(list1, list2_transformed, &scale_r2); STATUS("R2(I) = %5.4f %% (scale=%5.2e)\n", R2*100.0, scale_r2); - R1i = stat_r1_i(ref1, ref2_transformed, icommon, &scale_r1i); + R1i = stat_r1_i(list1, list2_transformed, &scale_r1i); STATUS("R1(I) = %5.4f %% (scale=%5.2e)\n", R1i*100.0, scale_r1i); - Rdiff = stat_rdiff_ignore(ref1, ref2_transformed, icommon, &scale_rdig); - STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (ignoring negative intensities)\n", + Rdiff = stat_rdiff_ignore(list1, list2_transformed, &scale_rdig); + STATUS("Rint(F) = %5.4f %% (scale=%5.2e)" + " (ignoring negative intensities)\n", Rdiff*100.0, scale_rdig); - Rdiff = stat_rdiff_zero(ref1, ref2_transformed, icommon, &scale); - STATUS("Rint(F) = %5.4f %% (scale=%5.2e) (zeroing negative intensities)\n", + Rdiff = stat_rdiff_zero(list1, list2_transformed, &scale); + STATUS("Rint(F) = %5.4f %% (scale=%5.2e)" + " (zeroing negative intensities)\n", Rdiff*100.0, scale); - Rdiff = stat_rdiff_intensity(ref1, ref2_transformed, icommon, - &scale_rintint); + Rdiff = stat_rdiff_intensity(list1, list2_transformed, &scale_rintint); STATUS("Rint(I) = %5.4f %% (scale=%5.2e)\n", Rdiff*100.0, scale_rintint); - pearson = stat_pearson_i(ref1, ref2_transformed, icommon); + pearson = stat_pearson_i(list1, list2_transformed); STATUS("Pearson r(I) = %5.4f\n", pearson); - pearson = stat_pearson_f_ignore(ref1, ref2_transformed, icommon); + pearson = stat_pearson_f_ignore(list1, list2_transformed); STATUS("Pearson r(F) = %5.4f (ignoring negative intensities)\n", pearson); - pearson = stat_pearson_f_zero(ref1, ref2_transformed, icommon); + pearson = stat_pearson_f_zero(list1, list2_transformed); STATUS("Pearson r(F) = %5.4f (zeroing negative intensities)\n", pearson); if ( config_shells ) { - plot_shells(ref1, ref2_transformed, icommon, scale_r1fi, + plot_shells(list1, list2_transformed, scale_r1fi, cell, sym, rmin_fix, rmax_fix); } - if ( outfile != NULL ) { - - write_reflections(outfile, icommon, out, NULL, - NULL, NULL, cell); - STATUS("Sigma(I) values in output file are not meaningful.\n"); - + if ( ratiofile != NULL ) { + write_reflist(ratiofile, ratio, cell); } return 0; diff --git a/src/get_hkl.c b/src/get_hkl.c index 1b7e4406..026198ab 100644 --- a/src/get_hkl.c +++ b/src/get_hkl.c @@ -22,7 +22,7 @@ #include <getopt.h> #include "utils.h" -#include "reflections.h" +#include "reflist-utils.h" #include "symmetry.h" #include "beam-parameters.h" @@ -45,7 +45,8 @@ static void show_help(const char *s) "To calculate Poisson samples accurately, you must also give:\n" " -b, --beam=<file> Get beam parameters from file.\n" "\n" -"You can artificially 'twin' the reflections, or expand them out:\n" +"You can artificially 'twin' the reflections, or expand them out. You can also" +" do both, in which case the 'twinning' will be done first:\n" " -w, --twin=<sym> Generate twinned data according to the given\n" " point group.\n" " -e, --expand=<sym> Expand reflections to this point group.\n" @@ -59,99 +60,123 @@ static void show_help(const char *s) "\n" "Don't forget to specify the output filename:\n" " -o, --output=<filename> Output filename (default: stdout).\n" +" -p, --pdb=<filename> Use unit cell parameters from this PDB file to\n" +" generate resolution values in the output file.\n" ); } /* Apply Poisson noise to all reflections */ -static void poisson_reflections(double *ref, ReflItemList *items, - double adu_per_photon) +static void poisson_reflections(RefList *list, double adu_per_photon) { - int i; - const int n = num_items(items); + Reflection *refl; + RefListIterator *iter; - for ( i=0; i<n; i++ ) { + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; - double val; - int c; + double val, c; - it = get_item(items, i); + val = get_intensity(refl); - val = lookup_intensity(ref, it->h, it->k, it->l); c = adu_per_photon * poisson_noise(val/adu_per_photon); - set_intensity(ref, it->h, it->k, it->l, c); - - progress_bar(i, n-1, "Simulating noise"); + set_int(refl, c); } } /* Apply 10% uniform noise to all reflections */ -static void noise_reflections(double *ref, ReflItemList *items) +static void noise_reflections(RefList *list) { - int i; - const int n = num_items(items); - - for ( i=0; i<n; i++ ) { + Reflection *refl; + RefListIterator *iter; - struct refl_item *it; - double val; - double r; + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - it = get_item(items, i); + double val, r; - val = lookup_intensity(ref, it->h, it->k, it->l); + val = get_intensity(refl); r = (double)random()/RAND_MAX; val += 0.1 * val * r; - set_intensity(ref, it->h, it->k, it->l, val); - - progress_bar(i, n-1, "Simulating noise"); + set_int(refl, val); } } -static ReflItemList *twin_reflections(double *ref, ReflItemList *items, - const char *holo, const char *mero, - double *esds) +static RefList *template_reflections(RefList *list, RefList *template) { - int i; - ReflItemList *new; + Reflection *refl; + RefListIterator *iter; + RefList *out; - new = new_items(); + out = reflist_new(); + + for ( refl = first_refl(template, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + signed int h, k, l; + Reflection *new; + Reflection *old; + + get_indices(refl, &h, &k, &l); + + old = find_refl(list, h, k, l); + if ( old == NULL ) continue; + + new = add_refl(out, h, k, l); + copy_data(new, old); + + } + + return out; +} + + +static RefList *twin_reflections(RefList *in, + const char *holo, const char *mero) +{ + Reflection *refl; + RefListIterator *iter; + RefList *out; if ( num_general_equivs(holo) < num_general_equivs(mero) ) { ERROR("%s is not a subgroup of %s!\n", mero, holo); return NULL; } - for ( i=0; i<num_items(items); i++ ) { + out = reflist_new(); + + for ( refl = first_refl(in, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { double total, sigma; - struct refl_item *it; signed int h, k, l; int n, j; int skip; - it = get_item(items, i); + get_indices(refl, &h, &k, &l); /* There is a many-to-one correspondence between reflections * in the merohedral and holohedral groups. Do the calculation * only once for each reflection in the holohedral group, which * contains fewer reflections. */ - get_asymm(it->h, it->k, it->l, &h, &k, &l, holo); - if ( find_item(new, h, k, l) ) continue; - - n = num_equivs(h, k, l, holo); + get_asymm(h, k, l, &h, &k, &l, holo); + if ( find_refl(out, h, k, l) != NULL ) continue; total = 0.0; sigma = 0.0; skip = 0; + n = num_equivs(h, k, l, holo); for ( j=0; j<n; j++ ) { signed int he, ke, le; @@ -164,8 +189,8 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items, * equivalent which belongs to our definition of the * asymmetric unit cell, so check them all. */ - if ( !find_unique_equiv(items, he, ke, le, mero, - &hu, &ku, &lu) ) { + if ( !find_equiv_in_list(in, he, ke, le, mero, + &hu, &ku, &lu) ) { /* Don't have this reflection, so bail out */ ERROR("Twinning %i %i %i requires the %i %i %i " "reflection (or an equivalent in %s), " @@ -176,61 +201,57 @@ static ReflItemList *twin_reflections(double *ref, ReflItemList *items, break; } - total += lookup_intensity(ref, hu, ku, lu); - sigma += pow(lookup_sigma(esds, hu, ku, lu), 2.0); + total += get_intensity(refl); + sigma += pow(get_esd_intensity(refl), 2.0); } if ( !skip ) { - set_intensity(ref, h, k, l, total); - set_sigma(esds, h, k, l, sqrt(sigma)); - add_item(new, h, k, l); + Reflection *new = add_refl(out, h, k, l); + set_int(new, total); + set_esd_intensity(new, sqrt(sigma)); } } - return new; + return out; } -static ReflItemList *expand_reflections(double *ref, ReflItemList *items, - const char *target, const char *initial) +static RefList *expand_reflections(RefList *in, + const char *target, const char *initial) { - int i; - ReflItemList *new; - - new = new_items(); + Reflection *refl; + RefListIterator *iter; + RefList *out; if ( num_general_equivs(target) > num_general_equivs(initial) ) { ERROR("%s is not a subgroup of %s!\n", initial, target); return NULL; } - for ( i=0; i<num_items(items); i++ ) { + out = reflist_new(); + + for ( refl = first_refl(in, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; signed int h, k, l; - signed int hd, kd, ld; int n, j; double intensity; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; - - /* Actually we don't really care what the equivalent is, - * we just want to be sure that there is nly be one version of - * this reflection. */ - find_unique_equiv(items, h, k, l, initial, &hd, &kd, &ld); + get_indices(refl, &h, &k, &l); + intensity = get_intensity(refl); - /* Now find out how many reflections need to be filled in */ n = num_equivs(h, k, l, initial); - intensity = lookup_intensity(ref, h, k, l); + /* For each equivalent in the higher symmetry group */ for ( j=0; j<n; j++ ) { signed int he, ke, le; + Reflection *new; /* Get the equivalent */ get_equiv(h, k, l, &he, &ke, &le, initial, j); @@ -239,43 +260,33 @@ static ReflItemList *expand_reflections(double *ref, ReflItemList *items, get_asymm(he, ke, le, &he, &ke, &le, target); /* Make sure the intensity is in the right place */ - set_intensity(ref, he, ke, le, intensity); - - /* Add the reflection, but only once */ - if ( !find_item(new, he, ke, le) ) { - add_item(new, he, ke, le); - } + new = add_refl(out, he, ke, le); + copy_data(new, refl); } } - return new; + return out; } int main(int argc, char *argv[]) { int c; - double *input_ref; - double *phases; - double *esds; - char *template = NULL; int config_noise = 0; int config_poisson = 0; - int config_nophase = 0; int config_multi = 0; char *holo = NULL; char *mero = NULL; char *expand = NULL; + char *input_file = NULL; + char *template = NULL; char *output = NULL; - char *input = NULL; - char *filename = NULL; - ReflItemList *input_items; - ReflItemList *write_items; - UnitCell *cell = NULL; char *beamfile = NULL; struct beam_params *beam = NULL; + RefList *input; + UnitCell *cell = NULL; /* Long options */ const struct option longopts[] = { @@ -290,6 +301,7 @@ int main(int argc, char *argv[]) {"intensities", 1, NULL, 'i'}, {"multiplicity", 0, &config_multi, 1}, {"beam", 1, NULL, 'b'}, + {"pdb", 1, NULL, 'p'}, {0, 0, NULL, 0} }; @@ -311,7 +323,7 @@ int main(int argc, char *argv[]) break; case 'i' : - input = strdup(optarg); + input_file = strdup(optarg); break; case 'y' : @@ -330,6 +342,14 @@ int main(int argc, char *argv[]) beamfile = strdup(optarg); break; + case 'p' : + cell = load_cell_from_pdb(optarg); + if ( cell == NULL ) { + ERROR("Failed to get cell from '%s'\n", optarg); + return 1; + } + break; + case 0 : break; @@ -339,10 +359,6 @@ int main(int argc, char *argv[]) } - if ( filename == NULL ) { - filename = strdup("molecule.pdb"); - } - if ( (holo != NULL) && (expand != NULL) ) { ERROR("You cannot 'twin' and 'expand' at the same time.\n"); ERROR("Decide which one you want to do first.\n"); @@ -358,19 +374,14 @@ int main(int argc, char *argv[]) } } - cell = load_cell_from_pdb(filename); - if ( !config_nophase ) { - phases = new_list_phase(); - } else { - phases = NULL; + if ( cell == NULL ) { + ERROR("You need to give a PDB file with the unit cell.\n"); + return 1; } - esds = new_list_sigma(); - input_ref = new_list_intensity(); - input_items = read_reflections(input, input_ref, phases, - NULL, esds); + input = read_reflections(input_file); free(input); - if ( check_symmetry(input_items, mero) ) { + if ( check_list_symmetry(input, mero) ) { ERROR("The input reflection list does not appear to" " have symmetry %s\n", mero); return 1; @@ -378,8 +389,7 @@ int main(int argc, char *argv[]) if ( config_poisson ) { if ( beam != NULL ) { - poisson_reflections(input_ref, input_items, - beam->adu_per_photon); + poisson_reflections(input, beam->adu_per_photon); } else { ERROR("You must give a beam parameters file in order" " to calculate Poisson noise.\n"); @@ -387,72 +397,69 @@ int main(int argc, char *argv[]) } } - if ( config_noise ) noise_reflections(input_ref, input_items); + if ( config_noise ) noise_reflections(input); if ( holo != NULL ) { - ReflItemList *new; + RefList *new; STATUS("Twinning from %s into %s\n", mero, holo); - new = twin_reflections(input_ref, input_items, - holo, mero, esds); - delete_items(input_items); - input_items = new; + new = twin_reflections(input, holo, mero); + + /* Replace old with new */ + reflist_free(input); + input = new; + + /* The symmetry of the list has changed */ + free(mero); + mero = holo; } if ( expand != NULL ) { - ReflItemList *new; + RefList *new; STATUS("Expanding from %s into %s\n", mero, expand); - new = expand_reflections(input_ref, input_items, expand, mero); - delete_items(input_items); - input_items = new; + new = expand_reflections(input, expand, mero); + + /* Replace old with new */ + reflist_free(input); + input = new; } if ( config_multi ) { - int i; + Reflection *refl; + RefListIterator *iter; - for ( i=0; i<num_items(input_items); i++ ) { + for ( refl = first_refl(input, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; double inty; + signed int h, k, l; + + get_indices(refl, &h, &k, &l); + inty = get_intensity(refl); - it = get_item(input_items, i); - inty = lookup_intensity(input_ref, it->h, it->k, it->l); - inty *= num_equivs(it->h, it->k, it->l, mero); - set_intensity(input_ref, it->h, it->k, it->l, inty); - STATUS("%i %i %i %i\n", it->h, it->k, it->l, - num_equivs(it->h, it->k, it->l, mero)); + inty *= (double)num_equivs(h, k, l, mero); + set_int(refl, inty); } } 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, NULL); - write_items = intersection_items(input_items, template_items); - delete_items(template_items); - - } else { - - /* Write out all reflections */ - write_items = new_items(); - /* (quick way of copying a list) */ - union_items(write_items, input_items); + RefList *t = read_reflections(template); + RefList *new = template_reflections(input, t); + reflist_free(input); + input = new; } - write_reflections(output, write_items, input_ref, esds, phases, - NULL, cell); + write_reflist(output, input, cell); - delete_items(input_items); - delete_items(write_items); + reflist_free(input); return 0; } diff --git a/src/indexamajig.c b/src/indexamajig.c index e1e0e8fd..c0757edd 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -31,11 +31,11 @@ #include "peaks.h" #include "detector.h" #include "filters.h" -#include "reflections.h" #include "thread-pool.h" #include "beam-parameters.h" #include "geometry.h" #include "stream.h" +#include "reflist-utils.h" enum { diff --git a/src/partialator.c b/src/partialator.c index 88926582..553a216a 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -26,7 +26,6 @@ #include "utils.h" #include "hdf5-file.h" #include "symmetry.h" -#include "reflections.h" #include "stream.h" #include "geometry.h" #include "peaks.h" @@ -35,6 +34,7 @@ #include "post-refinement.h" #include "hrs-scaling.h" #include "reflist.h" +#include "reflist-utils.h" static void show_help(const char *s) @@ -418,7 +418,7 @@ int main(int argc, char *argv[]) } /* Output results */ - write_reflections(outfile, obs, I_full, NULL, NULL, cts, NULL); + //write_reflist(outfile, obs, I_full, NULL, NULL, cts, NULL); /* Clean up */ for ( i=0; i<n_total_patterns; i++ ) { diff --git a/src/pattern_sim.c b/src/pattern_sim.c index 6d3c07ab..ba760e32 100644 --- a/src/pattern_sim.c +++ b/src/pattern_sim.c @@ -29,9 +29,10 @@ #include "hdf5-file.h" #include "detector.h" #include "peaks.h" -#include "reflections.h" #include "beam-parameters.h" #include "symmetry.h" +#include "reflist.h" +#include "reflist-utils.h" static void show_help(const char *s) @@ -412,34 +413,29 @@ int main(int argc, char *argv[]) } else { - int i; - ReflItemList *items; + RefList *reflections; + + reflections = read_reflections(intfile); + free(intfile); if ( grad == GRADIENT_PHASED ) { - phases = new_list_phase(); + phases = phases_from_list(reflections); } else { phases = NULL; } - intensities = new_list_intensity(); - phases = new_list_phase(); - flags = new_list_flag(); - items = read_reflections(intfile, intensities, phases, - NULL, NULL); - free(intfile); - - for ( i=0; i<num_items(items); i++ ) { - struct refl_item *it = get_item(items, i); - set_flag(flags, it->h, it->k, it->l, 1); - } + intensities = intensities_from_list(reflections); + phases = phases_from_list(reflections); + flags = flags_from_list(reflections); /* Check that the intensities have the correct symmetry */ - if ( check_symmetry(items, sym) ) { + if ( check_list_symmetry(reflections, sym) ) { ERROR("The input reflection list does not appear to" " have symmetry %s\n", sym); return 1; } - delete_items(items); + reflist_free(reflections); + } image.det = get_detector_geometry(geometry); diff --git a/src/povray.c b/src/povray.c index 808d3e22..8d0ec048 100644 --- a/src/povray.c +++ b/src/povray.c @@ -23,13 +23,13 @@ #include "utils.h" #include "symmetry.h" #include "render_hkl.h" +#include "povray.h" #define MAX_PROC (256) -int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts, - ReflItemList *items, unsigned int nproc, +int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc, const char *sym, int wght, double boost) { FILE *fh; @@ -39,6 +39,8 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts, pid_t pids[MAX_PROC]; float max; int i; + Reflection *refl; + RefListIterator *iter; if ( (nproc > MAX_PROC) || (nproc < 1) ) { ERROR("Number of processes must be a number between 1 and %i\n", @@ -165,27 +167,29 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts, fprintf(fh, "}\n"); max = 0.0; - for ( i=0; i<num_items(items); i++ ) { + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; float val; + signed int h, k, l; - it = get_item(items, i); + get_indices(refl, &h, &k, &l); switch ( wght ) { case WGHT_I : - val = lookup_intensity(ref, it->h, it->k, it->l); + val = get_intensity(refl); break; case WGHT_SQRTI : - val = lookup_intensity(ref, it->h, it->k, it->l); + val = get_intensity(refl); val = (val>0.0) ? sqrt(val) : 0.0; break; case WGHT_COUNTS : - val = lookup_count(counts, it->h, it->k, it->l); - val /= (float)num_equivs(it->h, it->k, it->l, sym); + val = get_redundancy(refl); + val /= (float)num_equivs(h, k, l, sym); break; case WGHT_RAWCOUNTS : - val = lookup_count(counts, it->h, it->k, it->l); + val = get_redundancy(refl); break; default : ERROR("Invalid weighting.\n"); @@ -197,30 +201,31 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts, } max /= boost; - for ( i=0; i<num_items(items); i++ ) { + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { - struct refl_item *it; - float radius; + signed int h, k, l;float radius; int s; float val, p, r, g, b, trans; int j; - it = get_item(items, i); + get_indices(refl, &h, &k, &l); switch ( wght ) { case WGHT_I : - val = lookup_intensity(ref, it->h, it->k, it->l); + val = get_intensity(refl); break; case WGHT_SQRTI : - val = lookup_intensity(ref, it->h, it->k, it->l); + val = get_intensity(refl); val = (val>0.0) ? sqrt(val) : 0.0; break; case WGHT_COUNTS : - val = lookup_count(counts, it->h, it->k, it->l); - val /= (float)num_equivs(it->h, it->k, it->l, sym); + val = get_redundancy(refl); + val /= (float)num_equivs(h, k, l, sym); break; case WGHT_RAWCOUNTS : - val = lookup_count(counts, it->h, it->k, it->l); + val = get_redundancy(refl); break; default : ERROR("Invalid weighting.\n"); @@ -269,12 +274,12 @@ int povray_render_animation(UnitCell *cell, double *ref, unsigned int *counts, trans = 1.0-(val/max); /* For each equivalent */ - for ( j=0; j<num_equivs(it->h, it->k, it->l, sym); j++ ) { + for ( j=0; j<num_equivs(h, k, l, sym); j++ ) { signed int he, ke, le; float x, y, z; - get_equiv(it->h, it->k, it->l, &he, &ke, &le, sym, j); + get_equiv(h, k, l, &he, &ke, &le, sym, j); x = asx*he + bsx*ke + csx*le; y = asy*he + bsy*ke + csy*le; diff --git a/src/povray.h b/src/povray.h index 69bf657e..0b0b3446 100644 --- a/src/povray.h +++ b/src/povray.h @@ -3,7 +3,7 @@ * * Invoke POV-ray * - * (c) 2006-2010 Thomas White <taw@physics.org> + * (c) 2006-2011 Thomas White <taw@physics.org> * * Part of CrystFEL - crystallography with a FEL * @@ -16,8 +16,10 @@ #ifndef POVRAY_H #define POVRAY_H -extern int povray_render_animation(UnitCell *cell, double *ref, - unsigned int *counts, ReflItemList *items, +#include "reflist.h" +#include "cell.h" + +extern int povray_render_animation(UnitCell *cell, RefList *list, unsigned int nproc, const char *sym, int wght, double boost); diff --git a/src/process_hkl.c b/src/process_hkl.c index 6ff40159..f5e705bd 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -468,11 +468,11 @@ int main(int argc, char *argv[]) output = strdup("processed.hkl"); } - if ( pdb == NULL ) { - pdb = strdup("molecule.pdb"); - } - cell = load_cell_from_pdb(pdb); + if ( cell == NULL ) { + ERROR("Failed to load cell from '%s'\n", pdb); + return 1; + } free(pdb); if ( sym == NULL ) sym = strdup("1"); @@ -537,7 +537,6 @@ int main(int argc, char *argv[]) write_reflist(output, model, cell); - cell_free(cell); fclose(fh); free(sym); reflist_free(model); diff --git a/src/reflections.c b/src/reflections.c deleted file mode 100644 index 4743aeee..00000000 --- a/src/reflections.c +++ /dev/null @@ -1,230 +0,0 @@ -/* - * reflections.c - * - * Utilities for handling reflections - * - * (c) 2006-2010 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - - -#include <stdlib.h> -#include <math.h> -#include <stdio.h> -#include <complex.h> -#include <string.h> - -#include "utils.h" -#include "cell.h" -#include "reflections.h" -#include "beam-parameters.h" - - -double *poisson_esds(double *intensities, ReflItemList *items, - double adu_per_photon) -{ - int i; - double *esds = new_list_sigma(); - - for ( i=0; i<num_items(items); i++ ) { - - struct refl_item *it; - signed int h, k, l; - double intensity, sigma; - - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; - - intensity = lookup_intensity(intensities, h, k, l); - - if ( intensity > 0.0 ) { - sigma = adu_per_photon * sqrt(intensity/adu_per_photon); - } else { - sigma = 0.0; - } - - set_sigma(esds, h, k, l, sigma); - - } - - return esds; -} - - -void write_reflections(const char *filename, ReflItemList *items, - double *intensities, double *esds, double *phases, - unsigned int *counts, UnitCell *cell) -{ - FILE *fh; - int i; - - if ( filename == NULL ) { - fh = stdout; - } else { - fh = fopen(filename, "w"); - } - - if ( fh == NULL ) { - ERROR("Couldn't open output file '%s'.\n", filename); - return; - } - - /* Write spacings and angle if zone axis pattern */ - fprintf(fh, " h k l I phase sigma(I) " - " 1/d(nm^-1) counts\n"); - - for ( i=0; i<num_items(items); i++ ) { - - struct refl_item *it; - signed int h, k, l; - - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; - - int N; - double intensity, s, sigma; - char ph[32]; - - if ( counts != NULL ) { - N = lookup_count(counts, h, k, l); - } else { - N = 1; - } - - intensity = lookup_intensity(intensities, h, k, l); - - if ( phases != NULL ) { - double p; - p = lookup_phase(phases, h, k, l); - snprintf(ph, 31, "%8.6f", p); - } else { - strncpy(ph, " -", 31); - } - - if ( cell != NULL ) { - s = 2.0*resolution(cell, h, k, l); - } else { - s = 0.0; - } - - if ( esds != NULL ) { - sigma = lookup_sigma(esds, h, k, l); - } else { - sigma = 0.0; - } - - /* h, k, l, I, sigma(I), s */ - fprintf(fh, "%3i %3i %3i %10.2f %s %10.2f %10.2f %7i\n", - h, k, l, intensity, ph, sigma, s/1.0e9, N); - - } - fclose(fh); -} - - -/* 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 *esds) -{ - FILE *fh; - char *rval; - ReflItemList *items; - - fh = fopen(filename, "r"); - if ( fh == NULL ) { - ERROR("Failed to open input file '%s'\n", filename); - return NULL; - } - - items = new_items(); - - do { - - char line[1024]; - signed int h, k, l; - float intensity, ph, res, sigma; - char phs[1024]; - int r; - int cts; - - rval = fgets(line, 1023, fh); - if ( rval == NULL ) continue; - r = sscanf(line, "%i %i %i %f %s %f %f %i", - &h, &k, &l, &intensity, phs, &sigma, &res, &cts); - if ( r >= 8 ) { - /* Woohoo */ - } else if ( r >= 7 ) { - /* No "counts", that's fine.. */ - cts = 1; - } else if ( r >= 6 ) { - /* No resolution. Didn't want it anyway. */ - res = 0.0; - } else if ( r >= 5 ) { - /* No sigma. It's OK today, but one - * day I'll get you... */ - sigma = 0.0; - } else if ( r >= 4 ) { - /* No phase. Better not need it.. */ - if ( phases != NULL ) { - ERROR("Need phases and none were specified!\n"); - abort(); - } - } else { - /* You lose. */ - continue; - } - - 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); - } - if ( esds != NULL ) { - set_sigma(esds, h, k, l, sigma); - } - - } while ( rval != NULL ); - - fclose(fh); - - return items; -} - - -double *ideal_intensities(double complex *sfac) -{ - double *ref; - signed int h, k, l; - - ref = new_list_intensity(); - - /* Generate ideal reflections from complex structure factors */ - for ( h=-INDMAX; h<=INDMAX; h++ ) { - for ( k=-INDMAX; k<=INDMAX; k++ ) { - for ( l=-INDMAX; l<=INDMAX; l++ ) { - double complex F = lookup_sfac(sfac, h, k, l); - double intensity = pow(cabs(F), 2.0); - set_intensity(ref, h, k, l, intensity); - } - } - } - - return ref; -} diff --git a/src/reflections.h b/src/reflections.h deleted file mode 100644 index 20060f21..00000000 --- a/src/reflections.h +++ /dev/null @@ -1,38 +0,0 @@ -/* - * reflections.h - * - * Utilities for handling reflections - * - * (c) 2006-2010 Thomas White <taw@physics.org> - * - * Part of CrystFEL - crystallography with a FEL - * - */ - -#ifdef HAVE_CONFIG_H -#include <config.h> -#endif - -#ifndef REFLECTIONS_H -#define REFLECTIONS_H - - -#include "cell.h" -#include "utils.h" - - -extern double *poisson_esds(double *intensities, ReflItemList *items, - double adu_per_photon); - -extern void write_reflections(const char *filename, ReflItemList *items, - double *intensities, double *esds, double *phases, - unsigned int *counts, UnitCell *cell); - -extern ReflItemList *read_reflections(const char *filename, - double *intensities, double *phases, - unsigned int *counts, double *esds); - -extern double *ideal_intensities(double complex *sfac); - - -#endif /* REFLECTIONS_H */ diff --git a/src/reflist-utils.c b/src/reflist-utils.c index 33d373b8..aac94fbf 100644 --- a/src/reflist-utils.c +++ b/src/reflist-utils.c @@ -15,6 +15,153 @@ #include "reflist.h" #include "cell.h" +#include "utils.h" +#include "reflist-utils.h" +#include "symmetry.h" + + +double *intensities_from_list(RefList *list) +{ + Reflection *refl; + RefListIterator *iter; + double *out = new_list_intensity(); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + signed int h, k, l; + double intensity = get_intensity(refl); + + get_indices(refl, &h, &k, &l); + + set_intensity(out, h, k, l, intensity); + + } + + return out; +} + + +double *phases_from_list(RefList *list) +{ + Reflection *refl; + RefListIterator *iter; + double *out = new_list_phase(); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + signed int h, k, l; + double phase = get_phase(refl); + + get_indices(refl, &h, &k, &l); + + set_phase(out, h, k, l, phase); + + } + + return out; + +} + + +unsigned char *flags_from_list(RefList *list) +{ + Reflection *refl; + RefListIterator *iter; + unsigned char *out = new_list_flag(); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + signed int h, k, l; + + get_indices(refl, &h, &k, &l); + + set_flag(out, h, k, l, 1); + + } + + return out; + +} + + +int check_list_symmetry(RefList *list, const char *sym) +{ + unsigned char *flags; + Reflection *refl; + RefListIterator *iter; + + flags = flags_from_list(list); + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + int j; + int found = 0; + signed int h, k, l; + + get_indices(refl, &h, &k, &l); + + for ( j=0; j<num_equivs(h, k, l, sym); j++ ) { + + signed int he, ke, le; + get_equiv(h, k, l, &he, &ke, &le, sym, j); + + if ( abs(he) > INDMAX ) continue; + if ( abs(le) > INDMAX ) continue; + if ( abs(ke) > INDMAX ) continue; + + found += lookup_flag(flags, he, ke, le); + + } + + if ( found > 1 ) { + free(flags); + return 1; /* Symmetry is wrong! */ + } + + } + + free(flags); + + return 0; +} + + +int find_equiv_in_list(RefList *list, signed int h, signed int k, + signed int l, const char *sym, signed int *hu, + signed int *ku, signed int *lu) +{ + int i; + int found = 0; + + for ( i=0; i<num_equivs(h, k, l, sym); i++ ) { + + signed int he, ke, le; + Reflection *f; + get_equiv(h, k, l, &he, &ke, &le, sym, i); + f = find_refl(list, he, ke, le); + + /* There must only be one equivalent. If there are more, it + * indicates that the user lied about the input symmetry. + * This situation should have been checked for earlier by + * calling check_symmetry() with 'items' and 'mero'. */ + + if ( (f != NULL) && !found ) { + *hu = he; *ku = ke; *lu = le; + return 1; + } + + } + + return 0; +} void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell) @@ -31,18 +178,26 @@ void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell) signed int h, k, l; double intensity, esd_i, s; + int red; double fs, ss; + char res[16]; 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! */ + esd_i = get_esd_intensity(refl); + red = get_redundancy(refl); + + if ( cell != NULL ) { + s = resolution(cell, h, k, l); + snprintf(res, 16, "%10.2f", s/1e9); + } else { + strcpy(res, " -"); + } - /* 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, + "%3i %3i %3i %10.2f %s %10.2f %s %7i %6.1f %6.1f\n", + h, k, l, intensity, " -", esd_i, res, red, fs, ss); } @@ -70,3 +225,85 @@ int write_reflist(const char *filename, RefList *list, UnitCell *cell) return 0; } + + +RefList *read_reflections_from_file(FILE *fh) +{ + char *rval = NULL; + int first = 1; + RefList *out; + + out = reflist_new(); + + do { + + char line[1024]; + signed int h, k, l; + float intensity, sigma, fs, ss; + char phs[1024]; + char ress[1024]; + int cts; + int r; + Reflection *refl; + + rval = fgets(line, 1023, fh); + if ( rval == NULL ) continue; + chomp(line); + + if ( strcmp(line, REFLECTION_END_MARKER) == 0 ) return out; + + r = sscanf(line, "%i %i %i %f %s %f %s %i %f %f", + &h, &k, &l, &intensity, phs, &sigma, ress, &cts, + &fs, &ss); + if ( (r != 10) && (!first) ) { + reflist_free(out); + return NULL; + } + + first = 0; + if ( r == 10 ) { + + double ph; + char *v; + + refl = add_refl(out, h, k, l); + set_int(refl, intensity); + set_detector_pos(refl, fs, ss, 0.0); + set_esd_intensity(refl, sigma); + + ph = strtod(phs, &v); + if ( v != NULL ) set_ph(refl, ph); + + /* The 1/d value is actually ignored. */ + + } + + } while ( rval != NULL ); + + /* Got read error of some kind before finding PEAK_LIST_END_MARKER */ + return NULL; +} + + +RefList *read_reflections(const char *filename) +{ + FILE *fh; + RefList *out; + + if ( filename == NULL ) { + fh = stdout; + } else { + fh = fopen(filename, "r"); + } + + if ( fh == NULL ) { + ERROR("Couldn't open input file '%s'.\n", filename); + return NULL; + } + + out = read_reflections_from_file(fh); + + fclose(fh); + + return out; +} diff --git a/src/reflist-utils.h b/src/reflist-utils.h index d5c1e655..b7ecd3ea 100644 --- a/src/reflist-utils.h +++ b/src/reflist-utils.h @@ -21,10 +21,24 @@ #include "cell.h" +#define REFLECTION_END_MARKER "End of reflections" + + extern void write_reflections_to_file(FILE *fh, RefList *list, UnitCell *cell); -extern int write_reflist(const char *filename, RefList *list, - UnitCell *cell); +extern int write_reflist(const char *filename, RefList *list, UnitCell *cell); + +extern RefList *read_reflections_from_file(FILE *fh); + +extern RefList *read_reflections(const char *filename); + +extern double *intensities_from_list(RefList *list); +extern double *phases_from_list(RefList *list); +extern unsigned char *flags_from_list(RefList *list); +extern int check_list_symmetry(RefList *list, const char *sym); +extern int find_equiv_in_list(RefList *list, signed int h, signed int k, + signed int l, const char *sym, signed int *hu, + signed int *ku, signed int *lu); #endif /* REFLIST_UTILS_H */ diff --git a/src/reflist.c b/src/reflist.c index 3181b418..a81c738e 100644 --- a/src/reflist.c +++ b/src/reflist.c @@ -46,6 +46,9 @@ struct _refldata { double intensity; double esd_i; + /* Phase */ + double phase; + /* Redundancy */ int redundancy; @@ -250,8 +253,20 @@ double get_esd_intensity(Reflection *refl) } +double get_phase(Reflection *refl) +{ + return refl->data.phase; +} + + /********************************** Setters ***********************************/ +void copy_data(Reflection *to, Reflection *from) +{ + memcpy(&to->data, &from->data, sizeof(struct _refldata)); +} + + void set_detector_pos(Reflection *refl, double exerr, double x, double y) { refl->data.excitation_error = exerr; @@ -301,6 +316,12 @@ void set_esd_intensity(Reflection *refl, double esd) } +void set_ph(Reflection *refl, double phase) +{ + refl->data.phase = phase; +} + + /********************************* Insertion **********************************/ static void insert_node(Reflection *head, Reflection *new) @@ -689,3 +710,9 @@ void optimise_reflist(RefList *list) recursive_depth(list->head->child[0])); } } + + +int num_reflections(RefList *list) +{ + return recursive_count(list->head->child[0]); +} diff --git a/src/reflist.h b/src/reflist.h index f6d01ecf..5b647102 100644 --- a/src/reflist.h +++ b/src/reflist.h @@ -44,8 +44,10 @@ 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); +extern double get_phase(Reflection *refl); /* Set */ +extern void copy_data(Reflection *to, Reflection *from); extern void set_detector_pos(Reflection *refl, double exerr, double x, double y); extern void set_partial(Reflection *refl, double r1, double r2, double p, @@ -55,6 +57,7 @@ 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); +extern void set_ph(Reflection *refl, double phase); /* Insertion */ extern Reflection *add_refl(RefList *list, INDICES); @@ -66,7 +69,8 @@ extern void delete_refl(Reflection *refl); extern Reflection *first_refl(RefList *list, RefListIterator **iterator); extern Reflection *next_refl(Reflection *refl, RefListIterator *iter); -/* Voodoo */ +/* Misc */ extern void optimise_reflist(RefList *list); +extern int num_reflections(RefList *list); #endif /* REFLIST_H */ diff --git a/src/render_hkl.c b/src/render_hkl.c index fc1c5121..aec6d714 100644 --- a/src/render_hkl.c +++ b/src/render_hkl.c @@ -26,11 +26,12 @@ #endif #include "utils.h" -#include "reflections.h" #include "povray.h" #include "symmetry.h" #include "render.h" #include "render_hkl.h" +#include "reflist.h" +#include "reflist-utils.h" static void show_help(const char *s) @@ -81,8 +82,7 @@ static void show_help(const char *s) static void draw_circles(signed int xh, signed int xk, signed int xl, signed int yh, signed int yk, signed int yl, signed int zh, signed int zk, signed int zl, - double *ref, unsigned int *counts, ReflItemList *items, - const char *sym, + RefList *list, const char *sym, cairo_t *dctx, int wght, double boost, int colscale, UnitCell *cell, double radius, double theta, double as, double bs, double cx, double cy, @@ -98,36 +98,36 @@ static void draw_circles(signed int xh, signed int xk, signed int xl, *max_res = 0.0; *max_ux = 0; *max_uy = 0; } - /* Loop across the two basis directions */ + /* Iterate over possible reflections in this zone */ for ( xi=-INDMAX; xi<INDMAX; xi++ ) { for ( yi=-INDMAX; yi<INDMAX; yi++ ) { double u, v, val, res; signed int h, k, l; - signed int he, ke, le; + Reflection *refl; h = xi*xh + yi*yh; k = xi*xk + yi*yk; l = xi*xl + yi*yl; /* Got this reflection? */ - if ( find_unique_equiv(items, h, k, l, sym, - &he, &ke, &le) == 0 ) continue; + refl = find_refl(list, h, k, l); + if ( refl == NULL ) continue; switch ( wght ) { case WGHT_I : - val = lookup_intensity(ref, he, ke, le); + val = get_intensity(refl); break; case WGHT_SQRTI : - val = lookup_intensity(ref, he, ke, le); + val = get_intensity(refl); val = (val>0.0) ? sqrt(val) : 0.0; break; case WGHT_COUNTS : - val = lookup_count(counts, he, ke, le); - val /= (float)num_equivs(he, ke, le, sym); + val = get_redundancy(refl); + val /= (float)num_equivs(h, k, l, sym); break; case WGHT_RAWCOUNTS : - val = lookup_count(counts, he, ke, le); + val = get_redundancy(refl); break; default : ERROR("Invalid weighting.\n"); @@ -228,8 +228,7 @@ static void render_overlined_indices(cairo_t *dctx, } -static void render_za(UnitCell *cell, ReflItemList *items, - double *ref, unsigned int *counts, +static void render_za(UnitCell *cell, RefList *list, double boost, const char *sym, int wght, int colscale, signed int xh, signed int xk, signed int xl, signed int yh, signed int yk, signed int yl, @@ -284,7 +283,7 @@ static void render_za(UnitCell *cell, ReflItemList *items, scale = 1.0; draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl, - ref, counts, items, sym, NULL, wght, boost, colscale, cell, + list, sym, NULL, wght, boost, colscale, cell, 0.0, theta, as, bs, 0.0, 0.0, scale, &max_ux, &max_uy, &max_val, &max_u, &max_v, &max_res); @@ -342,7 +341,7 @@ static void render_za(UnitCell *cell, ReflItemList *items, cy = 512.0 - 20.0; draw_circles(xh, xk, xl, yh, yk, yl, zh, zk, zl, - ref, counts, items, sym, dctx, wght, boost, colscale, cell, + list, sym, dctx, wght, boost, colscale, cell, max_r, theta, as, bs, cx, cy, scale, NULL, NULL, &max_val, NULL, NULL, NULL); @@ -462,8 +461,8 @@ int main(int argc, char *argv[]) { int c; UnitCell *cell; + RefList *list; char *infile; - double *ref; int config_povray = 0; int config_zoneaxis = 0; int config_sqrt = 0; @@ -477,7 +476,6 @@ int main(int argc, char *argv[]) int wght; int colscale; char *cscale = NULL; - unsigned int *cts; signed int dh=1, dk=0, dl=0; signed int rh=0, rk=1, rl=0; char *down = NULL; @@ -613,7 +611,8 @@ int main(int argc, char *argv[]) if ( config_zoneaxis ) { if ( (( down == NULL ) && ( right != NULL )) || (( down != NULL ) && ( right == NULL )) ) { - ERROR("Either specify both 'down' and 'right', or neither.\n"); + ERROR("Either specify both 'down' and 'right'," + " or neither.\n"); return 1; } if ( down != NULL ) { @@ -641,24 +640,22 @@ 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(); - ReflItemList *items = read_reflections(infile, ref, NULL, cts, NULL); - if ( items == NULL ) { - ERROR("Couldn't open file '%s'\n", infile); + list = read_reflections(infile); + if ( list == NULL ) { + ERROR("Couldn't read file '%s'\n", infile); return 1; } - if ( check_symmetry(items, sym) ) { + if ( check_list_symmetry(list, sym) ) { ERROR("The input reflection list does not appear to" " have symmetry %s\n", sym); return 1; } if ( config_povray ) { - r = povray_render_animation(cell, ref, cts, items, + r = povray_render_animation(cell, list, nproc, sym, wght, boost); } else if ( config_zoneaxis ) { - render_za(cell, items, ref, cts, boost, sym, wght, colscale, + render_za(cell, list, boost, sym, wght, colscale, rh, rk, rl, dh, dk, dl, outfile); } else { ERROR("Try again with either --povray or --zone-axis.\n"); @@ -666,7 +663,7 @@ int main(int argc, char *argv[]) free(pdb); free(sym); - delete_items(items); + reflist_free(list); if ( outfile != NULL ) free(outfile); return r; diff --git a/src/statistics.c b/src/statistics.c index d0d3ae85..77e7919f 100644 --- a/src/statistics.c +++ b/src/statistics.c @@ -25,9 +25,8 @@ struct r_params { - const double *ref1; - const double *ref2; - ReflItemList *items; /* Which reflections to use */ + RefList *list1; + RefList *list2; int fom; /* Which FoM to use (see the enum just below) */ }; @@ -43,27 +42,30 @@ enum { /* Return the least squares optimal scaling factor when comparing intensities. - * ref1,ref2 are the two intensity lists to compare. "items" is a ReflItemList + * list1,list2 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 stat_scale_intensity(RefList *list1, RefList *list2) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { double i1, i2; - struct refl_item *it; signed int h, k, l; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ - i1 = lookup_intensity(ref1, h, k, l); - i2 = lookup_intensity(ref2, h, k, l); + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); top += i1 * i2; bot += i2 * i2; @@ -77,30 +79,36 @@ double stat_scale_intensity(const double *ref1, const double *ref2, /* 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, + * list1,list2 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) +static double stat_scale_sqrti(RefList *list1, RefList *list2) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - double i1, i2, f1, f2; - struct refl_item *it; + double i1, i2; + double f1, f2; signed int h, k, l; + Reflection *refl2; + + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); - i1 = lookup_intensity(ref1, h, k, l); if ( i1 < 0.0 ) continue; f1 = sqrt(i1); - i2 = lookup_intensity(ref2, h, k, l); + if ( i2 < 0.0 ) continue; f2 = sqrt(i2); @@ -113,26 +121,33 @@ static double stat_scale_sqrti(const double *ref1, const double *ref2, } -static double internal_r1_ignorenegs(const double *ref1, const double *ref2, - ReflItemList *items, double scale) +static double internal_r1_ignorenegs(RefList *list1, RefList *list2, + double scale) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - double i1, f1, i2, f2; - struct refl_item *it; + double i1, i2; + double f1, f2; signed int h, k, l; + Reflection *refl2; + + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); - i1 = lookup_intensity(ref1, h, k, l); if ( i1 < 0.0 ) continue; f1 = sqrt(i1); - i2 = lookup_intensity(ref2, h, k, l); + if ( i2 < 0.0 ) continue; f2 = sqrt(i2); f2 *= scale; @@ -146,25 +161,32 @@ static double internal_r1_ignorenegs(const double *ref1, const double *ref2, } -static double internal_r1_negstozero(const double *ref1, const double *ref2, - ReflItemList *items, double scale) +static double internal_r1_negstozero(RefList *list1, RefList *list2, + double scale) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - double i1, f1, i2, f2; - struct refl_item *it; + double i1, i2; + double f1, f2; signed int h, k, l; + Reflection *refl2; + + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); - 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; f2 *= scale; @@ -177,24 +199,27 @@ static double internal_r1_negstozero(const double *ref1, const double *ref2, } -static double internal_r2(const double *ref1, const double *ref2, - ReflItemList *items, double scale) +static double internal_r2(RefList *list1, RefList *list2, double scale) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { double i1, i2; - struct refl_item *it; signed int h, k, l; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ - i1 = lookup_intensity(ref1, h, k, l); - i2 = scale * lookup_intensity(ref2, h, k, l); + i1 = get_intensity(refl1); + i2 = scale * get_intensity(refl2); top += pow(i1 - i2, 2.0); bot += pow(i1, 2.0); @@ -205,24 +230,27 @@ static double internal_r2(const double *ref1, const double *ref2, } -static double internal_r_i(const double *ref1, const double *ref2, - ReflItemList *items, double scale) +static double internal_r_i(RefList *list1, RefList *list2, double scale) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { double i1, i2; - struct refl_item *it; signed int h, k, l; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ - i1 = lookup_intensity(ref1, h, k, l); - i2 = scale * lookup_intensity(ref2, h, k, l); + i1 = get_intensity(refl1); + i2 = scale * get_intensity(refl2); top += fabs(i1-i2); bot += fabs(i1); @@ -233,24 +261,29 @@ static double internal_r_i(const double *ref1, const double *ref2, } -static double internal_rdiff_intensity(const double *ref1, const double *ref2, - ReflItemList *items, double scale) +static double internal_rdiff_intensity(RefList *list1, RefList *list2, + double scale) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { double i1, i2; - struct refl_item *it; signed int h, k, l; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ + + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); - i1 = lookup_intensity(ref1, h, k, l); - i2 = lookup_intensity(ref2, h, k, l); i2 *= scale; top += fabs(i1 - i2); @@ -262,25 +295,32 @@ static double internal_rdiff_intensity(const double *ref1, const double *ref2, } -static double internal_rdiff_negstozero(const double *ref1, const double *ref2, - ReflItemList *items, double scale) +static double internal_rdiff_negstozero(RefList *list1, RefList *list2, + double scale) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - double i1, i2, f1, f2; - struct refl_item *it; + double i1, i2; + double f1, f2; signed int h, k, l; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ + + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); - 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; f2 *= scale; @@ -293,26 +333,33 @@ static double internal_rdiff_negstozero(const double *ref1, const double *ref2, } -static double internal_rdiff_ignorenegs(const double *ref1, const double *ref2, - ReflItemList *items, double scale) +static double internal_rdiff_ignorenegs(RefList *list1, RefList *list2, + double scale) { double top = 0.0; double bot = 0.0; - int i; + Reflection *refl1; + RefListIterator *iter; - for ( i=0; i<num_items(items); i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - double i1, i2, f1, f2; - struct refl_item *it; + double i1, i2; + double f1, f2; signed int h, k, l; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ + + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); - i1 = lookup_intensity(ref1, h, k, l); if ( i1 < 0.0 ) continue; f1 = sqrt(i1); - i2 = lookup_intensity(ref2, h, k, l); + if ( i2 < 0.0 ) continue; f2 = sqrt(i2); f2 *= scale; @@ -332,26 +379,21 @@ static double calc_r(double scale, void *params) switch ( rp->fom) { case R_1_ZERO : - return internal_r1_negstozero(rp->ref1, rp->ref2, - rp->items, scale); + return internal_r1_negstozero(rp->list1, rp->list2, scale); case R_1_IGNORE : - return internal_r1_ignorenegs(rp->ref1, rp->ref2, - rp->items, scale); + return internal_r1_ignorenegs(rp->list1, rp->list2, scale); case R_2 : - return internal_r2(rp->ref1, rp->ref2, rp->items, scale); + return internal_r2(rp->list1, rp->list2, scale); case R_1_I : - return internal_r_i(rp->ref1, rp->ref2, rp->items, scale); + return internal_r_i(rp->list1, rp->list2, scale); case R_DIFF_ZERO : - return internal_rdiff_negstozero(rp->ref1, rp->ref2, - rp->items, scale); + return internal_rdiff_negstozero(rp->list1, rp->list2,scale); case R_DIFF_IGNORE : - return internal_rdiff_ignorenegs(rp->ref1, rp->ref2, - rp->items, scale); + return internal_rdiff_ignorenegs(rp->list1, rp->list2, scale); case R_DIFF_INTENSITY : - return internal_rdiff_intensity(rp->ref1, rp->ref2, - rp->items, scale); + return internal_rdiff_intensity(rp->list1, rp->list2, scale); } ERROR("No such FoM!\n"); @@ -359,8 +401,8 @@ static double calc_r(double scale, void *params) } -static double r_minimised(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep, int fom) +static double r_minimised(RefList *list1, RefList *list2, + double *scalep, int fom) { gsl_function F; gsl_min_fminimizer *s; @@ -369,9 +411,8 @@ static double r_minimised(const double *ref1, const double *ref2, struct r_params rp; int iter = 0; - rp.ref1 = ref1; - rp.ref2 = ref2; - rp.items = items; + rp.list1 = list1; + rp.list2 = list2; rp.fom = fom; F.function = &calc_r; @@ -385,12 +426,12 @@ static double r_minimised(const double *ref1, const double *ref2, case R_1_IGNORE : case R_DIFF_ZERO : case R_DIFF_IGNORE : - scale = stat_scale_sqrti(ref1, ref2, items); + scale = stat_scale_sqrti(list1, list2); break; case R_2 : case R_1_I : case R_DIFF_INTENSITY : - scale = stat_scale_intensity(ref1, ref2, items); + scale = stat_scale_intensity(list1, list2); break; } //STATUS("Initial scale factor estimate: %5.2e\n", scale); @@ -429,77 +470,74 @@ static double r_minimised(const double *ref1, const double *ref2, } -double stat_r1_ignore(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep) +double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep) { - return r_minimised(ref1, ref2, items, scalep, R_1_IGNORE); + return r_minimised(list1, list2, scalep, R_1_IGNORE); } -double stat_r1_zero(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep) +double stat_r1_zero(RefList *list1, RefList *list2, double *scalep) { - return r_minimised(ref1, ref2, items, scalep, R_1_ZERO); + return r_minimised(list1, list2, scalep, R_1_ZERO); } -double stat_r2(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep) +double stat_r2(RefList *list1, RefList *list2, double *scalep) { - return r_minimised(ref1, ref2, items, scalep, R_2); + return r_minimised(list1, list2, scalep, R_2); } -double stat_r1_i(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep) +double stat_r1_i(RefList *list1, RefList *list2, double *scalep) { - return r_minimised(ref1, ref2, items, scalep, R_1_I); + return r_minimised(list1, list2, scalep, R_1_I); } -double stat_rdiff_zero(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep) +double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep) { - return r_minimised(ref1, ref2, items, scalep, R_DIFF_ZERO); + return r_minimised(list1, list2, scalep, R_DIFF_ZERO); } -double stat_rdiff_ignore(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep) +double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep) { - return r_minimised(ref1, ref2, items, scalep, R_DIFF_IGNORE); + return r_minimised(list1, list2, scalep, R_DIFF_IGNORE); } -double stat_rdiff_intensity(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep) +double stat_rdiff_intensity(RefList *list1, RefList *list2, double *scalep) { - return r_minimised(ref1, ref2, items, scalep, R_DIFF_INTENSITY); + return r_minimised(list1, list2, scalep, R_DIFF_INTENSITY); } -double stat_pearson_i(const double *ref1, const double *ref2, - ReflItemList *items) + +double stat_pearson_i(RefList *list1, RefList *list2) { double *vec1, *vec2; - int i = 0; - int ni = num_items(items); + int ni = num_reflections(list1) + num_reflections(list2); double val; int nacc = 0; + Reflection *refl1; + RefListIterator *iter; vec1 = malloc(ni*sizeof(double)); vec2 = malloc(ni*sizeof(double)); - for ( i=0; i<ni; i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { double i1, i2; - struct refl_item *it; signed int h, k, l; + Reflection *refl2; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ - i1 = lookup_intensity(ref1, h, k, l); - i2 = lookup_intensity(ref2, h, k, l); + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); vec1[nacc] = i1; vec2[nacc] = i2; @@ -514,40 +552,42 @@ double stat_pearson_i(const double *ref1, const double *ref2, } -double stat_pearson_f_ignore(const double *ref1, const double *ref2, - ReflItemList *items) +double stat_pearson_f_ignore(RefList *list1, RefList *list2) { double *vec1, *vec2; - int i = 0; - int ni = num_items(items); + int ni = num_reflections(list1) + num_reflections(list2); double val; int nacc = 0; + Reflection *refl1; + RefListIterator *iter; vec1 = malloc(ni*sizeof(double)); vec2 = malloc(ni*sizeof(double)); - for ( i=0; i<ni; i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - double i1, i2, f1, f2; - struct refl_item *it; + double i1, i2; + double f1, f2; signed int h, k, l; + Reflection *refl2; int ig = 0; - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ + + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); - i1 = lookup_intensity(ref1, h, k, l); - if ( i1 < 0.0 ) { - ig=1; - } + if ( i1 < 0.0 ) ig = 1; f1 = sqrt(i1); - i2 = lookup_intensity(ref2, h, k, l); - if ( i2 < 0.0 ) { - ig=1; - } + + if ( i2 < 0.0 ) ig = 1; f2 = sqrt(i2); - if ( ig ) continue; + if ( ig ) continue; vec1[nacc] = f1; vec2[nacc] = f2; @@ -563,30 +603,35 @@ double stat_pearson_f_ignore(const double *ref1, const double *ref2, } -double stat_pearson_f_zero(const double *ref1, const double *ref2, - ReflItemList *items) +double stat_pearson_f_zero(RefList *list1, RefList *list2) { double *vec1, *vec2; - int i = 0; - int ni = num_items(items); + int ni = num_reflections(list1) + num_reflections(list2); double val; int nacc = 0; + Reflection *refl1; + RefListIterator *iter; vec1 = malloc(ni*sizeof(double)); vec2 = malloc(ni*sizeof(double)); - for ( i=0; i<ni; i++ ) { + for ( refl1 = first_refl(list1, &iter); + refl1 != NULL; + refl1 = next_refl(refl1, iter) ) { - double i1, i2, f1, f2; - struct refl_item *it; + double i1, i2; + double f1, f2; signed int h, k, l; + Reflection *refl2; + + get_indices(refl1, &h, &k, &l); + refl2 = find_refl(list2, h, k, l); + if ( refl2 == NULL ) continue; /* No common reflection */ - it = get_item(items, i); - h = it->h; k = it->k; l = it->l; + i1 = get_intensity(refl1); + i2 = get_intensity(refl2); - 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; vec1[nacc] = f1; diff --git a/src/statistics.h b/src/statistics.h index 18603357..34ef2406 100644 --- a/src/statistics.h +++ b/src/statistics.h @@ -17,35 +17,25 @@ #define STATISTICS_H -#include "utils.h" - -extern double stat_scale_intensity(const double *ref1, const double *ref2, - ReflItemList *items); - -extern double stat_r1_zero(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep); -extern double stat_r1_ignore(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_r1_i(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep); - -extern double stat_rdiff_zero(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep); -extern double stat_rdiff_ignore(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep); -extern double stat_rdiff_intensity(const double *ref1, const double *ref2, - ReflItemList *items, double *scalep); - -extern double stat_pearson_i(const double *ref1, const double *ref2, - ReflItemList *items); -extern double stat_pearson_f_zero(const double *ref1, const double *ref2, - ReflItemList *items); -extern double stat_pearson_f_ignore(const double *ref1, const double *ref2, - ReflItemList *items); +#include "reflist.h" + +extern double stat_scale_intensity(RefList *list1, RefList *list2); + +extern double stat_r1_zero(RefList *list1, RefList *list2, double *scalep); +extern double stat_r1_ignore(RefList *list1, RefList *list2, double *scalep); + +extern double stat_r2(RefList *list1, RefList *list2, double *scalep); + +extern double stat_r1_i(RefList *list1, RefList *list2, double *scalep); + +extern double stat_rdiff_zero(RefList *list1, RefList *list2, double *scalep); +extern double stat_rdiff_ignore(RefList *list1, RefList *list2, double *scalep); +extern double stat_rdiff_intensity(RefList *list1, RefList *list2, + double *scalep); + +extern double stat_pearson_i(RefList *list1, RefList *list2); +extern double stat_pearson_f_zero(RefList *list1, RefList *list2); +extern double stat_pearson_f_ignore(RefList *list1, RefList *list2); #endif /* STATISTICS_H */ diff --git a/src/stream.c b/src/stream.c index 545bed43..97c38b77 100644 --- a/src/stream.c +++ b/src/stream.c @@ -31,7 +31,8 @@ #define PEAK_LIST_START_MARKER "Peaks from peak search" #define PEAK_LIST_END_MARKER "End of peak list" #define REFLECTION_START_MARKER "Reflections measured after indexing" -#define REFLECTION_END_MARKER "End of reflections" +/* REFLECTION_END_MARKER is over in reflist-utils.h because it is also + * used to terminate a standalone list of reflections */ static void exclusive(const char *a, const char *b) { @@ -132,49 +133,6 @@ int count_patterns(FILE *fh) } -static int read_reflections(FILE *fh, struct image *image) -{ - char *rval = NULL; - int first = 1; - - image->reflections = reflist_new(); - - do { - - char line[1024]; - signed int h, k, l; - float intensity, sigma, res, fs, ss; - char phs[1024]; - int cts; - int r; - Reflection *refl; - - rval = fgets(line, 1023, fh); - if ( rval == NULL ) continue; - chomp(line); - - 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) && (!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 ); - - /* Got read error of some kind before finding PEAK_LIST_END_MARKER */ - return 1; -} - - static int read_peaks(FILE *fh, struct image *image) { char *rval = NULL; @@ -409,7 +367,8 @@ int read_chunk(FILE *fh, struct image *image) } if ( strcmp(line, REFLECTION_START_MARKER) == 0 ) { - if ( read_reflections(fh, image) ) { + image->reflections = read_reflections_from_file(fh); + if ( image->reflections == NULL ) { ERROR("Failed while reading reflections\n"); return 1; } |