1 files changed, 199 insertions, 123 deletions

diff --git a/src/im-sandbox.c b/src/im-sandbox.c
index 70839f21..aeeb023c 100644
--- a/src/im-sandbox.c
+++ b/src/im-sandbox.c
@@ -3,13 +3,13 @@
  *
  * Sandbox for indexing
  *
- * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY,
- *                  a research centre of the Helmholtz Association.
+ * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY,
+ *                       a research centre of the Helmholtz Association.
  * Copyright © 2012 Richard Kirian
  * Copyright © 2012 Lorenzo Galli
  *
  * Authors:
- *   2010-2012 Thomas White <taw@physics.org>
+ *   2010-2013 Thomas White <taw@physics.org>
  *   2011      Richard Kirian
  *   2012      Lorenzo Galli
  *   2012      Chunhong Yoon
@@ -47,6 +47,7 @@
 #include <sys/wait.h>
 #include <fcntl.h>
 #include <signal.h>
+#include <sys/stat.h>
 
 #ifdef HAVE_CLOCK_GETTIME
 #include <time.h>
@@ -73,14 +74,27 @@
 #define STATS_EVERY_N_SECONDS (5)
 
 
+/* Information about the indexing process for one pattern */
+struct pattern_args
+{
+	/* "Input" */
+	char *filename;
+
+	/* "Output" */
+	int n_crystals;
+};
+
+
 struct sandbox
 {
 	pthread_mutex_t lock;
 
-	int n_indexable;
 	int n_processed;
-	int n_indexable_last_stats;
+	int n_hadcrystals;
+	int n_crystals;
 	int n_processed_last_stats;
+	int n_hadcrystals_last_stats;
+	int n_crystals_last_stats;
 	int t_last_stats;
 
 	struct index_args *iargs;
@@ -91,6 +105,7 @@ struct sandbox
 	FILE **fhs;
 
 	int *running;
+	int *waiting;
 	FILE **result_fhs;
 	int *filename_pipes;
 	int *stream_pipe_read;
@@ -164,32 +179,30 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead,
 
 
 static void process_image(const struct index_args *iargs,
-                          struct pattern_args *pargs, FILE *ofh,
+                          struct pattern_args *pargs, Stream *st,
                           int cookie)
 {
 	float *data_for_measurement;
 	size_t data_size;
-	UnitCell *cell = iargs->cell;
-	int config_cmfilter = iargs->config_cmfilter;
-	int config_noisefilter = iargs->config_noisefilter;
-	int config_verbose = iargs->config_verbose;
-	IndexingMethod *indm = iargs->indm;
-	struct beam_params *beam = iargs->beam;
 	int check;
 	struct hdfile *hdfile;
 	struct image image;
+	int i;
+	char filename[1024];
+
+	/* Prefix to jump out of temporary folder */
+	snprintf(filename, 1023, "../../%s", pargs->filename);
 
 	image.features = NULL;
 	image.data = NULL;
 	image.flags = NULL;
-	image.indexed_cell = NULL;
-	image.det = copy_geom(iargs->det);
 	image.copyme = iargs->copyme;
-	image.reflections = NULL;
-	image.n_saturated = 0;
 	image.id = cookie;
-	image.filename = pargs->filename;
-	image.beam = beam;
+	image.filename = filename;
+	image.beam = iargs->beam;
+	image.det = iargs->det;
+	image.crystals = NULL;
+	image.n_crystals = 0;
 
 	hdfile = hdfile_open(image.filename);
 	if ( hdfile == NULL ) return;
@@ -231,53 +244,31 @@ static void process_image(const struct index_args *iargs,
 		      image.width, image.height,
 		      image.det->max_fs + 1, image.det->max_ss + 1);
 		hdfile_close(hdfile);
-		free_detector_geometry(image.det);
 		return;
 	}
 
-	if ( image.lambda < 0.0 ) {
-		if ( beam != NULL ) {
-			ERROR("Using nominal photon energy of %.2f eV\n",
-			      beam->photon_energy);
-			image.lambda = ph_en_to_lambda(
-			                          eV_to_J(beam->photon_energy));
-		} else {
-			ERROR("No wavelength in file, so you need to give "
-				"a beam parameters file with -b.\n");
-			hdfile_close(hdfile);
-			free_detector_geometry(image.det);
-			return;
-		}
-	}
-
-	if ( image.lambda > 1000 ) {
-		if ( beam != NULL ) {
-			ERROR("Nonsensical wavelength in HDF5."
-			      "Using nominal photon energy of %.2f eV\n",
-			      beam->photon_energy);
-			image.lambda = ph_en_to_lambda(
-			                          eV_to_J(beam->photon_energy));
-		} else {
-			ERROR("Nonsensical wavelength in file, so you need to "
-			      "give a beam parameters file with -b.\n");
-			hdfile_close(hdfile);
-			free_detector_geometry(image.det);
-			return;
-		}
-	}
-
 	fill_in_values(image.det, hdfile);
+	fill_in_beam_parameters(image.beam, hdfile);
 
-	if ( config_cmfilter ) {
-		filter_cm(&image);
+	image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy));
+
+	if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) {
+		/* Error message covers a silly value in the beam file or in
+		 * the HDF5 file. */
+		ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n",
+		      image.lambda, image.beam->photon_energy, image.filename);
+		hdfile_close(hdfile);
+		return;
 	}
 
+	if ( iargs->cmfilter ) filter_cm(&image);
+
 	/* Take snapshot of image after CM subtraction but before
 	 * the aggressive noise filter. */
 	data_size = image.width * image.height * sizeof(float);
 	data_for_measurement = malloc(data_size);
 
-	if ( config_noisefilter ) {
+	if ( iargs->noisefilter ) {
 		filter_noise(&image, data_for_measurement);
 	} else {
 		memcpy(data_for_measurement, image.data, data_size);
@@ -291,14 +282,18 @@ static void process_image(const struct index_args *iargs,
 			iargs->hdf5_peak_path)) {
 			ERROR("Failed to get peaks from HDF5 file.\n");
 		}
-		validate_peaks(&image, iargs->min_int_snr,
-		               iargs->ir_inn, iargs->ir_mid, iargs->ir_out);
+		if ( !iargs->no_revalidate ) {
+			validate_peaks(&image, iargs->min_int_snr,
+				       iargs->ir_inn, iargs->ir_mid,
+				       iargs->ir_out, iargs->use_saturated);
+		}
 		break;
 
 		case PEAK_ZAEF:
 		search_peaks(&image, iargs->threshold,
 		             iargs->min_gradient, iargs->min_snr,
-		             iargs->ir_inn, iargs->ir_mid, iargs->ir_out);
+		             iargs->ir_inn, iargs->ir_mid, iargs->ir_out,
+		             iargs->use_saturated);
 		break;
 
 	}
@@ -308,50 +303,63 @@ static void process_image(const struct index_args *iargs,
 	free(image.data);
 	image.data = data_for_measurement;
 
-	/* Calculate orientation matrix (by magic) */
-	image.div = beam->divergence;
-	image.bw = beam->bandwidth;
-	image.profile_radius = beam->profile_radius;
-
-	index_pattern(&image, cell, indm, iargs->cellr,
-	              config_verbose, iargs->ipriv,
-	              iargs->config_insane, iargs->tols);
-
-	if ( image.indexed_cell != NULL ) {
-		pargs->indexable = 1;
-		image.reflections = find_intersections(&image,
-				image.indexed_cell);
-		if (image.reflections != NULL) {
-			integrate_reflections(&image,
-			                      iargs->config_closer,
-			                      iargs->config_bgsub,
-			                      iargs->min_int_snr,
-			                      iargs->ir_inn,
-			                      iargs->ir_mid,
-			                      iargs->ir_out);
+	/* Index the pattern */
+	index_pattern(&image, iargs->indm, iargs->ipriv);
+
+	pargs->n_crystals = image.n_crystals;
+
+	/* Default beam parameters */
+	image.div = image.beam->divergence;
+	image.bw = image.beam->bandwidth;
+
+	/* Integrate each crystal's diffraction spots */
+	for ( i=0; i<image.n_crystals; i++ ) {
+
+		RefList *reflections;
+
+		/* Set default crystal parameter(s) */
+		crystal_set_profile_radius(image.crystals[i],
+		                           image.beam->profile_radius);
+
+		if ( iargs->integrate_found ) {
+			reflections = select_intersections(&image,
+			                                   image.crystals[i]);
+		} else {
+			reflections = find_intersections(&image,
+			                                 image.crystals[i]);
 		}
-	} else {
-		image.reflections = NULL;
-	}
 
-	write_chunk(ofh, &image, hdfile, iargs->stream_flags);
-	fprintf(ofh, "END\n");
-	fflush(ofh);
+		crystal_set_reflections(image.crystals[i], reflections);
+
+	}
 
-	/* Only free cell if found */
-	cell_free(image.indexed_cell);
+	/* Integrate all the crystals at once - need all the crystals so that
+	 * overlaps can be detected. */
+	integrate_reflections(&image, iargs->closer,
+	                              iargs->bgsub,
+	                              iargs->min_int_snr,
+	                              iargs->ir_inn,
+	                              iargs->ir_mid,
+	                              iargs->ir_out,
+	                              iargs->integrate_saturated,
+	                              iargs->res_cutoff);
+
+	write_chunk(st, &image, hdfile,
+	            iargs->stream_peaks, iargs->stream_refls);
+
+	for ( i=0; i<image.n_crystals; i++ ) {
+		crystal_free(image.crystals[i]);
+	}
 
-	reflist_free(image.reflections);
 	free(image.data);
 	if ( image.flags != NULL ) free(image.flags);
 	image_feature_list_free(image.features);
 	hdfile_close(hdfile);
-	free_detector_geometry(image.det);
 }
 
 
 static void run_work(const struct index_args *iargs,
-                     int filename_pipe, int results_pipe, FILE *ofh,
+                     int filename_pipe, int results_pipe, Stream *st,
                      int cookie)
 {
 	int allDone = 0;
@@ -397,12 +405,12 @@ static void run_work(const struct index_args *iargs,
 		} else {
 
 			pargs.filename = line;
-			pargs.indexable = 0;
+			pargs.n_crystals = 0;
 
-			process_image(iargs, &pargs, ofh, cookie);
+			process_image(iargs, &pargs, st, cookie);
 
 			/* Request another image */
-			c = sprintf(buf, "%i\n", pargs.indexable);
+			c = sprintf(buf, "%i\n", pargs.n_crystals);
 			w = write(results_pipe, buf, c);
 			if ( w < 0 ) {
 				ERROR("write P0\n");
@@ -414,11 +422,13 @@ static void run_work(const struct index_args *iargs,
 
 	}
 
-	cleanup_indexing(iargs->ipriv);
+	write_line(st, "DONE");
+
+	cleanup_indexing(iargs->indm, iargs->ipriv);
 	free(iargs->indm);
 	free(iargs->ipriv);
 	free_detector_geometry(iargs->det);
-	free(iargs->beam);
+	free_beam_parameters(iargs->beam);
 	free(iargs->element);
 	free(iargs->hdf5_peak_path);
 	free_copy_hdf5_field_list(iargs->copyme);
@@ -479,15 +489,26 @@ static int pump_chunk(FILE *fh, FILE *ofh)
 
 		}
 
-		if ( strcmp(line, "END\n") == 0 ) {
+		if ( strcmp(line, "FLUSH\n") == 0 ) {
 			chunk_finished = 1;
-		} else {
+			continue;
+		}
+
+		if ( strcmp(line, "DONE\n") == 0 ) {
+			return 1;
+		}
+
+		fprintf(ofh, "%s", line);
+
+		if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) {
+			chunk_finished = 1;
+		}
+		if ( strcmp(line, CHUNK_START_MARKER"\n") == 0 ) {
 			chunk_started = 1;
-			fprintf(ofh, "%s", line);
 		}
 
-	} while ( !chunk_finished );
 
+	} while ( !chunk_finished );
 	return 0;
 }
 
@@ -541,10 +562,13 @@ static void *run_reader(void *sbv)
 
 			if ( !sb->running[i] ) continue;
 
-			if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) continue;
+			if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) {
+				continue;
+			}
 
 			if ( pump_chunk(sb->fhs[i], sb->ofh) ) {
 				sb->running[i] = 0;
+				sb->waiting[i] = 0;
 			}
 
 		}
@@ -561,7 +585,37 @@ static void *run_reader(void *sbv)
 }
 
 
-static void start_worker_process(struct sandbox *sb, int slot)
+static int create_temporary_folder(signed int n)
+{
+	int r;
+	char tmp[64];
+
+	if ( n < 0 ) {
+		snprintf(tmp, 63, "indexamajig.%i", getpid());
+	} else {
+		snprintf(tmp, 63, "worker.%i", n);
+	}
+
+	r = mkdir(tmp, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
+	if ( r ) {
+		ERROR("Failed to create temporary folder: %s\n",
+		      strerror(errno));
+		return 1;
+	}
+
+	r = chdir(tmp);
+	if ( r ) {
+		ERROR("Failed to chdir to temporary folder: %s\n",
+		      strerror(errno));
+		return 1;
+	}
+
+	return 0;
+}
+
+
+static void start_worker_process(struct sandbox *sb, int slot,
+                                 int argc, char *argv[])
 {
 	pid_t p;
 	int filename_pipe[2];
@@ -585,7 +639,7 @@ static void start_worker_process(struct sandbox *sb, int slot)
 
 	if ( p == 0 ) {
 
-		FILE *sfh;
+		Stream *st;
 		int j;
 		struct sigaction sa;
 		int r;
@@ -600,6 +654,8 @@ static void start_worker_process(struct sandbox *sb, int slot)
 			return;
 		}
 
+		create_temporary_folder(slot);
+
 		/* Free resources which will not be needed by worker */
 		for ( j=0; j<sb->n_proc; j++ ) {
 			if ( (j != slot) && (sb->running[j]) ) {
@@ -622,10 +678,12 @@ static void start_worker_process(struct sandbox *sb, int slot)
 		close(filename_pipe[1]);
 		close(result_pipe[0]);
 
-		sfh = fdopen(sb->stream_pipe_write[slot], "w");
+		st = open_stream_fd_for_write(sb->stream_pipe_write[slot]);
+		write_command(st, argc, argv);
+		write_line(st, "FLUSH");
 		run_work(sb->iargs, filename_pipe[0], result_pipe[1],
-		         sfh, slot);
-		fclose(sfh);
+		         st, slot);
+		close_stream(st);
 
 		//close(filename_pipe[0]);
 		close(result_pipe[1]);
@@ -671,6 +729,7 @@ static void handle_zombie(struct sandbox *sb)
 		int status, p;
 
 		if ( !sb->running[i] ) continue;
+		if ( sb->waiting[i] ) continue;
 
 		p = waitpid(sb->pids[i], &status, WNOHANG);
 
@@ -682,6 +741,7 @@ static void handle_zombie(struct sandbox *sb)
 		if ( p == sb->pids[i] ) {
 
 			sb->running[i] = 0;
+			sb->waiting[i] = 1;
 
 			if ( WIFEXITED(status) ) {
 				continue;
@@ -693,7 +753,7 @@ static void handle_zombie(struct sandbox *sb)
 				STATUS("Last filename was: %s\n",
 				       sb->last_filename[i]);
 				sb->n_processed++;
-				start_worker_process(sb, i);
+				start_worker_process(sb, i, 0, NULL);
 			}
 
 		}
@@ -705,7 +765,7 @@ static void handle_zombie(struct sandbox *sb)
 
 void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
                     int config_basename, FILE *fh, char *use_this_one_instead,
-                    FILE *ofh)
+                    FILE *ofh, int argc, char *argv[])
 {
 	int i;
 	int allDone;
@@ -720,17 +780,19 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		return;
 	}
 
-	sb->n_indexable = 0;
 	sb->n_processed = 0;
-	sb->n_indexable_last_stats = 0;
+	sb->n_hadcrystals = 0;
+	sb->n_crystals = 0;
 	sb->n_processed_last_stats = 0;
+	sb->n_hadcrystals_last_stats = 0;
+	sb->n_crystals_last_stats = 0;
 	sb->t_last_stats = get_monotonic_seconds();
 	sb->n_proc = n_proc;
-	sb->ofh = ofh;
 	sb->iargs = iargs;
 
 	pthread_mutex_init(&sb->lock, NULL);
 
+	sb->ofh = ofh;
 	sb->stream_pipe_read = calloc(n_proc, sizeof(int));
 	sb->stream_pipe_write = calloc(n_proc, sizeof(int));
 	if ( sb->stream_pipe_read == NULL ) {
@@ -761,6 +823,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 	sb->result_fhs = calloc(n_proc, sizeof(FILE *));
 	sb->pids = calloc(n_proc, sizeof(pid_t));
 	sb->running = calloc(n_proc, sizeof(int));
+	sb->waiting = calloc(n_proc, sizeof(int));
 	sb->fhs = calloc(sb->n_proc, sizeof(FILE *));
 	if ( sb->filename_pipes == NULL ) {
 		ERROR("Couldn't allocate memory for pipes.\n");
@@ -778,6 +841,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		ERROR("Couldn't allocate memory for process flags.\n");
 		return;
 	}
+	if ( sb->waiting == NULL ) {
+		ERROR("Couldn't allocate memory for process flags.\n");
+		return;
+	}
 
 	sb->last_filename = calloc(n_proc, sizeof(char *));
 	if ( sb->last_filename == NULL ) {
@@ -810,10 +877,12 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		return;
 	}
 
+	if ( create_temporary_folder(-1) ) return;
+
 	/* Fork the right number of times */
 	lock_sandbox(sb);
 	for ( i=0; i<n_proc; i++ ) {
-		start_worker_process(sb, i);
+		start_worker_process(sb, i, argc, argv);
 	}
 	unlock_sandbox(sb);
 
@@ -852,10 +921,9 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 
 		r = select(fdmax+1, &fds, NULL, NULL, &tv);
 		if ( r == -1 ) {
-			if ( errno != EINTR ) {
-				ERROR("select() failed: %s\n", strerror(errno));
-			}
-			continue;
+			if ( errno == EINTR ) continue;
+			ERROR("select() failed: %s\n", strerror(errno));
+			break;
 		}
 
 		if ( r == 0 ) continue; /* No progress this time.  Try again */
@@ -903,8 +971,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 					ERROR("Invalid result '%s'\n", results);
 				}
 			} else {
-				sb->n_indexable += atoi(results);
+
+				int nc = atoi(results);
+				sb->n_crystals += nc;
+				if ( nc > 0 ) {
+					sb->n_hadcrystals++;
+				}
 				sb->n_processed++;
+
 			}
 
 			/* Send next filename */
@@ -937,14 +1011,16 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 		tNow = get_monotonic_seconds();
 		if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) {
 
-			STATUS("%i out of %i indexed so far,"
-			       " %i out of %i since the last message.\n",
-			       sb->n_indexable, sb->n_processed,
-			       sb->n_indexable - sb->n_indexable_last_stats,
+			STATUS("%4i indexable out of %4i processed, "
+			       "%4i crystals so far. "
+			       "%4i images processed since the last message.\n",
+			       sb->n_hadcrystals, sb->n_processed,
+			       sb->n_crystals,
 			       sb->n_processed - sb->n_processed_last_stats);
 
-			sb->n_indexable_last_stats = sb->n_indexable;
 			sb->n_processed_last_stats = sb->n_processed;
+			sb->n_hadcrystals_last_stats = sb->n_hadcrystals;
+			sb->n_crystals_last_stats = sb->n_crystals;
 			sb->t_last_stats = tNow;
 
 		}
@@ -982,10 +1058,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix,
 	free(sb->result_fhs);
 	free(sb->pids);
 	free(sb->running);
+	free(sb->waiting);
 
-	if ( ofh != stdout ) fclose(ofh);
-
-	STATUS("There were %i images, of which %i could be indexed.\n",
-	       sb->n_processed, sb->n_indexable);
+	STATUS("Final:"
+	       " %i images processed, %i had crystals, %i crystals overall.\n",
+	       sb->n_processed, sb->n_hadcrystals, sb->n_crystals);
 
 }