facetron: Placeholders for HRS scaling

3 files changed, 43 insertions, 111 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index 5599da98..6031e441 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -64,7 +64,7 @@ calibrate_detector_LDADD = @LIBS@
 
 facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
                    image.c geometry.c reflections.c stream.c thread-pool.c \
-                   beam-parameters.c symmetry.c post-refinement.c
+                   beam-parameters.c symmetry.c post-refinement.c hrs-scaling.c
 facetron_LDADD = @LIBS@
 
 cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \
@@ -87,4 +87,5 @@ EXTRA_DIST = cell.h hdf5-file.h image.h  utils.h diffraction.h detector.h \
              render.h hdfsee.h dirax.h peaks.h index.h  filters.h \
              diffraction-gpu.h cl-utils.h symmetry.h \
              povray.h index-priv.h geometry.h templates.h render_hkl.h \
-             stream.h thread-pool.h beam-parameters.h post-refinement.h
+             stream.h thread-pool.h beam-parameters.h post-refinement.h \
+             hrs-scaling.h
diff --git a/src/Makefile.in b/src/Makefile.in
index fbdcd53c..a427a4ef 100644
--- a/src/Makefile.in
+++ b/src/Makefile.in
@@ -85,7 +85,7 @@ am_facetron_OBJECTS = facetron.$(OBJEXT) cell.$(OBJEXT) \
 	peaks.$(OBJEXT) image.$(OBJEXT) geometry.$(OBJEXT) \
 	reflections.$(OBJEXT) stream.$(OBJEXT) thread-pool.$(OBJEXT) \
 	beam-parameters.$(OBJEXT) symmetry.$(OBJEXT) \
-	post-refinement.$(OBJEXT)
+	post-refinement.$(OBJEXT) hrs-scaling.$(OBJEXT)
 facetron_OBJECTS = $(am_facetron_OBJECTS)
 facetron_DEPENDENCIES =
 am_get_hkl_OBJECTS = get_hkl.$(OBJEXT) sfac.$(OBJEXT) cell.$(OBJEXT) \
@@ -325,7 +325,7 @@ calibrate_detector_SOURCES = calibrate_detector.c utils.c hdf5-file.c image.c \
 calibrate_detector_LDADD = @LIBS@
 facetron_SOURCES = facetron.c cell.c hdf5-file.c utils.c detector.c peaks.c \
                    image.c geometry.c reflections.c stream.c thread-pool.c \
-                   beam-parameters.c symmetry.c post-refinement.c
+                   beam-parameters.c symmetry.c post-refinement.c hrs-scaling.c
 
 facetron_LDADD = @LIBS@
 cubeit_SOURCES = cubeit.c cell.c hdf5-file.c utils.c detector.c render.c \
@@ -346,7 +346,8 @@ EXTRA_DIST = cell.h hdf5-file.h image.h  utils.h diffraction.h detector.h \
              render.h hdfsee.h dirax.h peaks.h index.h  filters.h \
              diffraction-gpu.h cl-utils.h symmetry.h \
              povray.h index-priv.h geometry.h templates.h render_hkl.h \
-             stream.h thread-pool.h beam-parameters.h post-refinement.h
+             stream.h thread-pool.h beam-parameters.h post-refinement.h \
+             hrs-scaling.h
 
 all: all-am
 
@@ -487,6 +488,7 @@ distclean-compile:
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/get_hkl.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/hdf5-file.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/hdfsee.Po@am__quote@
+@AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/hrs-scaling.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/image.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/index.Po@am__quote@
 @AMDEP_TRUE@@am__include@ @am__quote@./$(DEPDIR)/indexamajig.Po@am__quote@
diff --git a/src/facetron.c b/src/facetron.c
index d14141c9..38d701e6 100644
--- a/src/facetron.c
+++ b/src/facetron.c
@@ -33,6 +33,7 @@
 #include "thread-pool.h"
 #include "beam-parameters.h"
 #include "post-refinement.h"
+#include "hrs-scaling.h"
 
 
 /* Maximum number of iterations of NLSq to do for each image per macrocycle. */
@@ -134,26 +135,39 @@ static void refine_image(int mytask, void *tasks)
 }
 
 
-struct integrate_args
+static void refine_all(struct image *images, int n_total_patterns,
+                       struct detector *det, const char *sym,
+                       ReflItemList *obs, double *i_full, int nthreads,
+                       FILE *graph, FILE *pgraph)
 {
-	const char *sym;
-	ReflItemList *obs;
-	double *i_full;
-	unsigned int *cts;
-	pthread_mutex_t *list_lock;
-	struct image *image;
-};
+	struct refine_args *tasks;
+	int i;
 
+	tasks = malloc(n_total_patterns * sizeof(struct refine_args));
+	for ( i=0; i<n_total_patterns; i++ ) {
+
+		tasks[i].sym = sym;
+		tasks[i].obs = obs;
+		tasks[i].i_full = i_full;
+		tasks[i].image = &images[i];
+		tasks[i].graph = graph;
+		tasks[i].pgraph = pgraph;
+
+	}
 
-static void integrate_image(int mytask, void *tasks)
+	run_thread_range(n_total_patterns, nthreads, "Refining",
+	                 refine_image, tasks);
+
+	free(tasks);
+}
+
+
+static void integrate_image(struct image *image)
 {
-	struct integrate_args *all_args = tasks;
-	struct integrate_args *pargs = &all_args[mytask];
 	struct cpeak *spots;
 	int j, n;
 	struct hdfile *hdfile;
-	struct image *image = pargs->image;
-	double nominal_photon_energy = pargs->image->beam->photon_energy;
+	double nominal_photon_energy = image->beam->photon_energy;
 
 	hdfile = hdfile_open(image->filename);
 	if ( hdfile == NULL ) {
@@ -165,7 +179,7 @@ static void integrate_image(int mytask, void *tasks)
 		return;
 	}
 
-	if ( hdf5_read(hdfile, pargs->image, 0, nominal_photon_energy) ) {
+	if ( hdf5_read(hdfile, image, 0, nominal_photon_energy) ) {
 		ERROR("Couldn't read '%s'\n", image->filename);
 		hdfile_close(hdfile);
 		return;
@@ -178,10 +192,8 @@ static void integrate_image(int mytask, void *tasks)
 	for ( j=0; j<n; j++ ) {
 
 		signed int h, k, l;
-		signed int ha, ka, la;
 		float i_partial;
 		float xc, yc;
-		float i_full_est;
 
 		h = spots[j].h;
 		k = spots[j].k;
@@ -198,26 +210,13 @@ static void integrate_image(int mytask, void *tasks)
 		                    &xc, &yc, &i_partial, NULL, NULL, 1, 1) ) {
 			continue;
 		}
-		i_partial *= image->osf;
-
-		i_full_est = i_partial / spots[j].p;
-
-		get_asymm(h, k, l, &ha, &ka, &la, pargs->sym);
-
-		pthread_mutex_lock(pargs->list_lock);
-		integrate_intensity(pargs->i_full, ha, ka, la, i_full_est);
-		integrate_count(pargs->cts, ha, ka, la, 1);
-		if ( !find_item(pargs->obs, ha, ka, la) ) {
-			add_item(pargs->obs, ha, ka, la);
-		}
-		pthread_mutex_unlock(pargs->list_lock);
 
 	}
 
 	free(image->data);
 	if ( image->flags != NULL ) free(image->flags);
 	hdfile_close(hdfile);
-	free(spots);
+	image->cpeaks = spots;
 
 	/* Muppet proofing */
 	image->data = NULL;
@@ -225,78 +224,6 @@ static void integrate_image(int mytask, void *tasks)
 }
 
 
-static void refine_all(struct image *images, int n_total_patterns,
-                       struct detector *det, const char *sym,
-                       ReflItemList *obs, double *i_full, int nthreads,
-                       FILE *graph, FILE *pgraph)
-{
-	struct refine_args *tasks;
-	int i;
-
-	tasks = malloc(n_total_patterns * sizeof(struct refine_args));
-	for ( i=0; i<n_total_patterns; i++ ) {
-
-		tasks[i].sym = sym;
-		tasks[i].obs = obs;
-		tasks[i].i_full = i_full;
-		tasks[i].image = &images[i];
-		tasks[i].graph = graph;
-		tasks[i].pgraph = pgraph;
-
-	}
-
-	run_thread_range(n_total_patterns, nthreads, "Refining",
-	                 refine_image, tasks);
-
-	free(tasks);
-}
-
-
-static void estimate_full(struct image *images, int n_total_patterns,
-                          struct detector *det, const char *sym,
-                          ReflItemList *obs, double *i_full, unsigned int *cts,
-                          int nthreads)
-{
-	int i;
-	struct integrate_args *tasks;
-	pthread_mutex_t list_lock = PTHREAD_MUTEX_INITIALIZER;
-
-	clear_items(obs);
-	memset(i_full, 0, LIST_SIZE*sizeof(double));
-	memset(cts, 0, LIST_SIZE*sizeof(unsigned int));
-
-	tasks = malloc(n_total_patterns * sizeof(struct integrate_args));
-	for ( i=0; i<n_total_patterns; i++ ) {
-
-		tasks[i].sym = sym;
-		tasks[i].obs = obs;
-		tasks[i].i_full = i_full;
-		tasks[i].cts = cts;
-		tasks[i].list_lock = &list_lock;
-		tasks[i].image = &images[i];
-
-	}
-
-	run_thread_range(n_total_patterns, nthreads, "Integrating",
-	                 integrate_image, tasks);
-
-	free(tasks);
-
-	/* Divide the totals to get the means */
-	for ( i=0; i<num_items(obs); i++ ) {
-
-		struct refl_item *it;
-		double total;
-
-		it = get_item(obs, i);
-		total = lookup_intensity(i_full, it->h, it->k, it->l);
-		total /= lookup_count(cts, it->h, it->k, it->l);
-		set_intensity(i_full, it->h, it->k, it->l, total);
-
-	}
-}
-
-
 int main(int argc, char *argv[])
 {
 	int c;
@@ -485,6 +412,9 @@ int main(int argc, char *argv[])
 		images[i].data = NULL;
 		images[i].flags = NULL;
 
+		/* Get reflections from this image */
+		integrate_image(&images[i]);
+
 		free(filename);
 
 		progress_bar(i, n_total_patterns-1, "Loading pattern data");
@@ -496,8 +426,8 @@ int main(int argc, char *argv[])
 	cts = new_list_count();
 
 	/* Make initial estimates */
-	estimate_full(images, n_total_patterns, det, sym, obs, i_full, cts,
-	              nthreads);
+	STATUS("Performing initial scaling.\n");
+	scale_intensities(images, n_total_patterns, sym);
 
 	/* Iterate */
 	for ( i=0; i<n_iter; i++ ) {
@@ -527,8 +457,7 @@ int main(int argc, char *argv[])
 		           nthreads, fhg, fhp);
 
 		/* Re-estimate all the full intensities */
-		estimate_full(images, n_total_patterns, det, sym, obs, i_full,
-		              cts, nthreads);
+		scale_intensities(images, n_total_patterns, sym);
 
 		fclose(fhg);
 		fclose(fhp);