aboutsummaryrefslogtreecommitdiff
path: root/src/indexamajig.c
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2012-07-23 15:00:33 +0200
committerThomas White <taw@physics.org>2012-07-23 15:00:33 +0200
commit9964cdeafc1338f5841270bc5392ad8e3c8eb33e (patch)
tree0fe784c1f27967034690b32bc67cb4426d79aae3 /src/indexamajig.c
parente63626f4207fa58d3ea70b68f2862630d281598d (diff)
parent6eb751c2230226723884cc7bc473d65b91064b81 (diff)
Merge branch 'tom/speed'
Diffstat (limited to 'src/indexamajig.c')
-rw-r--r--src/indexamajig.c540
1 files changed, 56 insertions, 484 deletions
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 1ff392e9..294bce1d 100644
--- a/src/indexamajig.c
+++ b/src/indexamajig.c
@@ -12,6 +12,7 @@
* 2010-2012 Thomas White <taw@physics.org>
* 2011 Richard Kirian
* 2012 Lorenzo Galli
+ * 2012 Chunhong Yoon
*
* This file is part of CrystFEL.
*
@@ -43,7 +44,6 @@
#include <getopt.h>
#include <hdf5.h>
#include <gsl/gsl_errno.h>
-#include <pthread.h>
#ifdef HAVE_CLOCK_GETTIME
#include <time.h>
@@ -63,81 +63,7 @@
#include "stream.h"
#include "reflist-utils.h"
-
-/* Write statistics at APPROXIMATELY this interval */
-#define STATS_EVERY_N_SECONDS (5)
-
-
-enum {
- PEAK_ZAEF,
- PEAK_HDF5,
-};
-
-
-/* Information about the indexing process which is common to all patterns */
-struct static_index_args
-{
- UnitCell *cell;
- int config_cmfilter;
- int config_noisefilter;
- int config_verbose;
- int stream_flags; /* What goes into the output? */
- int config_satcorr;
- int config_closer;
- int config_insane;
- int config_bgsub;
- float threshold;
- float min_gradient;
- float min_snr;
- double min_int_snr;
- struct detector *det;
- IndexingMethod *indm;
- IndexingPrivate **ipriv;
- int peaks; /* Peak detection method */
- int cellr;
- float tols[4];
- struct beam_params *beam;
- const char *element;
- const char *hdf5_peak_path;
- double ir_inn;
- double ir_mid;
- double ir_out;
-
- /* Output stream */
- pthread_mutex_t *output_mutex; /* Protects the output stream */
- FILE *ofh;
- const struct copy_hdf5_field *copyme;
-};
-
-
-/* Information about the indexing process for one pattern */
-struct index_args
-{
- /* "Input" */
- char *filename;
- struct static_index_args static_args;
-
- /* "Output" */
- int indexable;
-};
-
-
-/* Information needed to choose the next task and dispatch it */
-struct queue_args
-{
- FILE *fh;
- char *prefix;
- int config_basename;
- struct static_index_args static_args;
-
- char *use_this_one_instead;
-
- int n_indexable;
- int n_processed;
- int n_indexable_last_stats;
- int n_processed_last_stats;
- int t_last_stats;
-};
+#include "im-sandbox.h"
static void show_help(const char *s)
@@ -234,294 +160,10 @@ static void show_help(const char *s)
" least 10%% of the located peaks.\n"
" --no-bg-sub Don't subtract local background estimates from\n"
" integrated intensities.\n"
-"\n"
-"\nYou can tune the CPU affinities for enhanced performance on NUMA machines:\n"
-"\n"
-" --cpus=<n> Specify number of CPUs. This is NOT the same as\n"
-" giving the number of analyses to run in parallel.\n"
-" --cpugroup=<n> Batch threads in groups of this size.\n"
-" --cpuoffset=<n> Start using CPUs at this group number.\n"
);
}
-static void process_image(void *pp, int cookie)
-{
- struct index_args *pargs = pp;
- struct hdfile *hdfile;
- struct image image;
- float *data_for_measurement;
- size_t data_size;
- char *filename = pargs->filename;
- UnitCell *cell = pargs->static_args.cell;
- int config_cmfilter = pargs->static_args.config_cmfilter;
- int config_noisefilter = pargs->static_args.config_noisefilter;
- int config_verbose = pargs->static_args.config_verbose;
- IndexingMethod *indm = pargs->static_args.indm;
- struct beam_params *beam = pargs->static_args.beam;
-
- image.features = NULL;
- image.data = NULL;
- image.flags = NULL;
- image.indexed_cell = NULL;
- image.reflections = NULL;
- image.id = cookie;
- image.filename = filename;
- image.det = copy_geom(pargs->static_args.det);
- image.copyme = pargs->static_args.copyme;
- image.beam = beam;
-
- pargs->indexable = 0;
-
- hdfile = hdfile_open(filename);
- if ( hdfile == NULL ) return;
-
- if ( pargs->static_args.element != NULL ) {
-
- int r;
- r = hdfile_set_image(hdfile, pargs->static_args.element);
- if ( r ) {
- ERROR("Couldn't select path '%s'\n",
- pargs->static_args.element);
- hdfile_close(hdfile);
- return;
- }
-
- } else {
-
- int r;
- r = hdfile_set_first_image(hdfile, "/");
- if ( r ) {
- ERROR("Couldn't select first path\n");
- hdfile_close(hdfile);
- return;
- }
-
- }
-
- hdf5_read(hdfile, &image, pargs->static_args.config_satcorr);
-
- if ( (image.width != image.det->max_fs+1)
- || (image.height != image.det->max_ss+1) )
- {
- ERROR("Image size doesn't match geometry size"
- " - rejecting image.\n");
- ERROR("Image size: %i,%i. Geometry size: %i,%i\n",
- 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;
- }
- }
- fill_in_values(image.det, hdfile);
-
- if ( config_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 ) {
- filter_noise(&image, data_for_measurement);
- } else {
- memcpy(data_for_measurement, image.data, data_size);
- }
-
- switch ( pargs->static_args.peaks )
- {
- case PEAK_HDF5 :
- /* Get peaks from HDF5 */
- if ( get_peaks(&image, hdfile,
- pargs->static_args.hdf5_peak_path) )
- {
- ERROR("Failed to get peaks from HDF5 file.\n");
- }
- break;
-
- case PEAK_ZAEF :
- search_peaks(&image, pargs->static_args.threshold,
- pargs->static_args.min_gradient,
- pargs->static_args.min_snr,
- pargs->static_args.ir_inn,
- pargs->static_args.ir_mid,
- pargs->static_args.ir_out);
- break;
- }
-
- /* Get rid of noise-filtered version at this point
- * - it was strictly for the purposes of peak detection. */
- free(image.data);
- image.data = data_for_measurement;
-
- /* Calculate orientation matrix (by magic) */
- image.div = beam->divergence;
- image.bw = beam->bandwidth;
- image.profile_radius = 0.0001e9;
- index_pattern(&image, cell, indm, pargs->static_args.cellr,
- config_verbose, pargs->static_args.ipriv,
- pargs->static_args.config_insane,
- pargs->static_args.tols);
-
- if ( image.indexed_cell != NULL ) {
-
- pargs->indexable = 1;
-
- image.reflections = find_intersections(&image,
- image.indexed_cell);
-
- if ( image.reflections != NULL ) {
-
- integrate_reflections(&image,
- pargs->static_args.config_closer,
- pargs->static_args.config_bgsub,
- pargs->static_args.min_int_snr,
- pargs->static_args.ir_inn,
- pargs->static_args.ir_mid,
- pargs->static_args.ir_out);
-
- }
-
- } else {
-
- image.reflections = NULL;
-
- }
-
- pthread_mutex_lock(pargs->static_args.output_mutex);
- write_chunk(pargs->static_args.ofh, &image, hdfile,
- pargs->static_args.stream_flags);
- pthread_mutex_unlock(pargs->static_args.output_mutex);
-
- /* Only free cell if found */
- cell_free(image.indexed_cell);
-
- 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 *get_image(void *qp)
-{
- char *line;
- struct index_args *pargs;
- char *rval;
- struct queue_args *qargs = qp;
-
- /* Initialise new task arguments */
- pargs = malloc(sizeof(struct index_args));
- memcpy(&pargs->static_args, &qargs->static_args,
- sizeof(struct static_index_args));
-
- /* Get the next filename */
- if ( qargs->use_this_one_instead != NULL ) {
-
- line = qargs->use_this_one_instead;
- qargs->use_this_one_instead = NULL;
-
- } else {
-
- line = malloc(1024*sizeof(char));
- rval = fgets(line, 1023, qargs->fh);
- if ( rval == NULL ) {
- free(pargs);
- free(line);
- return NULL;
- }
- chomp(line);
-
- }
-
- if ( qargs->config_basename ) {
- char *tmp;
- tmp = safe_basename(line);
- free(line);
- line = tmp;
- }
-
- pargs->filename = malloc(strlen(qargs->prefix)+strlen(line)+1);
-
- snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line);
-
- free(line);
-
- return pargs;
-}
-
-
-#ifdef HAVE_CLOCK_GETTIME
-
-static time_t get_monotonic_seconds()
-{
- struct timespec tp;
- clock_gettime(CLOCK_MONOTONIC, &tp);
- return tp.tv_sec;
-}
-
-#else
-
-/* Fallback version of the above. The time according to gettimeofday() is not
- * monotonic, so measuring intervals based on it will screw up if there's a
- * timezone change (e.g. daylight savings) while the program is running. */
-static time_t get_monotonic_seconds()
-{
- struct timeval tp;
- gettimeofday(&tp, NULL);
- return tp.tv_sec;
-}
-
-#endif
-
-static void finalise_image(void *qp, void *pp)
-{
- struct queue_args *qargs = qp;
- struct index_args *pargs = pp;
- time_t monotonic_seconds;
-
- qargs->n_indexable += pargs->indexable;
- qargs->n_processed++;
-
- monotonic_seconds = get_monotonic_seconds();
- if ( monotonic_seconds >= qargs->t_last_stats+STATS_EVERY_N_SECONDS ) {
-
- STATUS("%i out of %i indexed so far,"
- " %i out of %i since the last message.\n",
- qargs->n_indexable, qargs->n_processed,
- qargs->n_indexable - qargs->n_indexable_last_stats,
- qargs->n_processed - qargs->n_processed_last_stats);
-
- qargs->n_processed_last_stats = qargs->n_processed;
- qargs->n_indexable_last_stats = qargs->n_indexable;
- qargs->t_last_stats = monotonic_seconds;
-
- }
-
- free(pargs->filename);
- free(pargs);
-}
-
-
static int parse_cell_reduction(const char *scellr, int *err,
int *reduction_needs_cell)
{
@@ -554,7 +196,6 @@ int main(int argc, char *argv[])
FILE *fh;
FILE *ofh;
char *rval = NULL;
- int n_images;
int config_noindex = 0;
int config_cmfilter = 0;
int config_noisefilter = 0;
@@ -585,19 +226,15 @@ int main(int argc, char *argv[])
float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */
int cellr;
int peaks;
- int nthreads = 1;
- pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER;
+ int n_proc = 1;
char *prepare_line;
char prepare_filename[1024];
- struct queue_args qargs;
+ char *use_this_one_instead;
+ struct index_args iargs;
struct beam_params *beam = NULL;
char *element = NULL;
double nominal_photon_energy;
int stream_flags = STREAM_INTEGRATED;
- int cpu_num = 0;
- int cpu_groupsize = 1;
- int cpu_offset = 0;
- char *endptr;
char *hdf5_peak_path = NULL;
struct copy_hdf5_field *copyme;
char *intrad = NULL;
@@ -684,7 +321,7 @@ int main(int argc, char *argv[])
break;
case 'j' :
- nthreads = atoi(optarg);
+ n_proc = atoi(optarg);
break;
case 'g' :
@@ -726,39 +363,10 @@ int main(int argc, char *argv[])
break;
case 6 :
- cpu_num = strtol(optarg, &endptr, 10);
- if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
- ERROR("Invalid number of CPUs ('%s')\n",
- optarg);
- return 1;
- }
- break;
-
case 7 :
- cpu_groupsize = strtol(optarg, &endptr, 10);
- if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
- ERROR("Invalid CPU group size ('%s')\n",
- optarg);
- return 1;
- }
- if ( cpu_groupsize < 1 ) {
- ERROR("CPU group size cannot be"
- " less than 1.\n");
- return 1;
- }
- break;
-
case 8 :
- cpu_offset = strtol(optarg, &endptr, 10);
- if ( !( (optarg[0] != '\0') && (endptr[0] == '\0') ) ) {
- ERROR("Invalid CPU offset ('%s')\n",
- optarg);
- return 1;
- }
- if ( cpu_offset < 0 ) {
- ERROR("CPU offset must be positive.\n");
- return 1;
- }
+ ERROR("The options --cpus, --cpugroup and --cpuoffset"
+ " are no longer used by indexamajig.\n");
break;
case 9 :
@@ -795,12 +403,6 @@ int main(int argc, char *argv[])
}
- if ( (cpu_num > 0) && (cpu_num % cpu_groupsize != 0) ) {
- ERROR("Number of CPUs must be divisible by"
- " the CPU group size.\n");
- return 1;
- }
-
if ( filename == NULL ) {
filename = strdup("-");
}
@@ -816,18 +418,15 @@ int main(int argc, char *argv[])
free(filename);
if ( outfile == NULL ) {
- outfile = strdup("-");
- }
- if ( strcmp(outfile, "-") == 0 ) {
ofh = stdout;
} else {
ofh = fopen(outfile, "w");
+ if ( ofh == NULL ) {
+ ERROR("Failed to open output file '%s'\n", outfile);
+ return 1;
+ }
+ free(outfile);
}
- if ( ofh == NULL ) {
- ERROR("Failed to open output file '%s'\n", outfile);
- return 1;
- }
- free(outfile);
if ( hdf5_peak_path == NULL ) {
hdf5_peak_path = strdup("/processing/hitfinder/peakinfo");
@@ -860,8 +459,8 @@ int main(int argc, char *argv[])
}
}
- if ( nthreads == 0 ) {
- ERROR("Invalid number of threads.\n");
+ if ( n_proc == 0 ) {
+ ERROR("Invalid number of processes.\n");
return 1;
}
@@ -929,6 +528,7 @@ int main(int argc, char *argv[])
ERROR("Invalid parameters for '--int-radius'\n");
return 1;
}
+ free(intrad);
} else {
STATUS("WARNING: You did not specify --int-radius.\n");
STATUS("WARNING: I will use the default values, which are"
@@ -967,22 +567,20 @@ int main(int argc, char *argv[])
} else {
STATUS("No beam parameters file was given, so I'm taking the"
" nominal photon energy to be 2 keV.\n");
+ ERROR("I'm also going to assume 1 ADU per photon, which is");
+ ERROR(" almost certainly wrong. Peak sigmas will be"
+ " incorrect.\n");
nominal_photon_energy = 2000.0;
}
- if ( beam == NULL ) {
- ERROR("Warning: no beam parameters file.\n");
- ERROR("I'm going to assume 1 ADU per photon, which is almost");
- ERROR(" certainly wrong. Peak sigmas will be incorrect.\n");
- }
-
/* Get first filename and use it to set up the indexing */
- prepare_line = malloc(1024*sizeof(char));
+ prepare_line = malloc(1024);
rval = fgets(prepare_line, 1023, fh);
if ( rval == NULL ) {
ERROR("Failed to get filename to prepare indexing.\n");
return 1;
}
+ use_this_one_instead = strdup(prepare_line);
chomp(prepare_line);
if ( config_basename ) {
char *tmp;
@@ -991,7 +589,7 @@ int main(int argc, char *argv[])
prepare_line = tmp;
}
snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
- qargs.use_this_one_instead = prepare_line;
+ free(prepare_line);
/* Prepare the indexer */
if ( indm != NULL ) {
@@ -1007,67 +605,41 @@ int main(int argc, char *argv[])
gsl_set_error_handler_off();
- qargs.static_args.cell = cell;
- qargs.static_args.config_cmfilter = config_cmfilter;
- qargs.static_args.config_noisefilter = config_noisefilter;
- qargs.static_args.config_verbose = config_verbose;
- qargs.static_args.config_satcorr = config_satcorr;
- qargs.static_args.config_closer = config_closer;
- qargs.static_args.config_insane = config_insane;
- qargs.static_args.config_bgsub = config_bgsub;
- qargs.static_args.cellr = cellr;
- qargs.static_args.tols[0] = tols[0];
- qargs.static_args.tols[1] = tols[1];
- qargs.static_args.tols[2] = tols[2];
- qargs.static_args.tols[3] = tols[3];
- qargs.static_args.threshold = threshold;
- qargs.static_args.min_gradient = min_gradient;
- qargs.static_args.min_snr = min_snr;
- qargs.static_args.min_int_snr = min_int_snr;
- qargs.static_args.det = det;
- qargs.static_args.indm = indm;
- qargs.static_args.ipriv = ipriv;
- qargs.static_args.peaks = peaks;
- qargs.static_args.output_mutex = &output_mutex;
- qargs.static_args.ofh = ofh;
- qargs.static_args.beam = beam;
- qargs.static_args.element = element;
- qargs.static_args.stream_flags = stream_flags;
- qargs.static_args.hdf5_peak_path = hdf5_peak_path;
- qargs.static_args.copyme = copyme;
- qargs.static_args.ir_inn = ir_inn;
- qargs.static_args.ir_mid = ir_mid;
- qargs.static_args.ir_out = ir_out;
-
- qargs.fh = fh;
- qargs.prefix = prefix;
- qargs.config_basename = config_basename;
- qargs.n_indexable = 0;
- qargs.n_processed = 0;
- qargs.n_indexable_last_stats = 0;
- qargs.n_processed_last_stats = 0;
- qargs.t_last_stats = get_monotonic_seconds();
-
- n_images = run_threads(nthreads, process_image, get_image,
- finalise_image, &qargs, 0,
- cpu_num, cpu_groupsize, cpu_offset);
-
- cleanup_indexing(ipriv);
-
- free(indm);
- free(ipriv);
+ /* Static worker args */
+ iargs.cell = cell;
+ iargs.config_cmfilter = config_cmfilter;
+ iargs.config_noisefilter = config_noisefilter;
+ iargs.config_verbose = config_verbose;
+ iargs.config_satcorr = config_satcorr;
+ iargs.config_closer = config_closer;
+ iargs.config_insane = config_insane;
+ iargs.config_bgsub = config_bgsub;
+ iargs.cellr = cellr;
+ iargs.tols[0] = tols[0];
+ iargs.tols[1] = tols[1];
+ iargs.tols[2] = tols[2];
+ iargs.tols[3] = tols[3];
+ iargs.threshold = threshold;
+ iargs.min_gradient = min_gradient;
+ iargs.min_snr = min_snr;
+ iargs.min_int_snr = min_int_snr;
+ iargs.det = det;
+ iargs.indm = indm;
+ iargs.ipriv = ipriv;
+ iargs.peaks = peaks;
+ iargs.beam = beam;
+ iargs.element = element;
+ iargs.stream_flags = stream_flags;
+ iargs.hdf5_peak_path = hdf5_peak_path;
+ iargs.copyme = copyme;
+ iargs.ir_inn = ir_inn;
+ iargs.ir_mid = ir_mid;
+ iargs.ir_out = ir_out;
+
+ create_sandbox(&iargs, n_proc, prefix, config_basename, fh,
+ use_this_one_instead, ofh);
+
free(prefix);
- free_detector_geometry(det);
- free(beam);
- free(element);
- free(hdf5_peak_path);
- free_copy_hdf5_field_list(copyme);
- cell_free(cell);
- if ( fh != stdin ) fclose(fh);
- if ( ofh != stdout ) fclose(ofh);
-
- STATUS("There were %i images, of which %i could be indexed.\n",
- n_images, qargs.n_indexable);
return 0;
}