diff options
-rw-r--r-- | .gitignore | 1 | ||||
-rw-r--r-- | doc/man/partial_sim.1 | 11 | ||||
-rw-r--r-- | libcrystfel/src/cell.c | 2 | ||||
-rw-r--r-- | libcrystfel/src/cell.h | 2 | ||||
-rw-r--r-- | libcrystfel/src/index.c | 2 | ||||
-rw-r--r-- | libcrystfel/src/index.h | 3 | ||||
-rw-r--r-- | libcrystfel/src/mosflm.c | 7 | ||||
-rw-r--r-- | libcrystfel/src/peaks.c | 7 | ||||
-rw-r--r-- | src/indexamajig.c | 738 | ||||
-rw-r--r-- | src/partial_sim.c | 4 | ||||
-rw-r--r-- | tests/symmetry_check.c | 5 |
11 files changed, 473 insertions, 309 deletions
@@ -46,3 +46,4 @@ build-aux/missing build-aux/config.guess build-aux/depcomp build-aux/ltmain.sh +/nbproject/private/
\ No newline at end of file diff --git a/doc/man/partial_sim.1 b/doc/man/partial_sim.1 index e41708ed..5b44d34f 100644 --- a/doc/man/partial_sim.1 +++ b/doc/man/partial_sim.1 @@ -71,17 +71,6 @@ When combined with with \fB-i\fR, specifies the symmetry of the input reflection .PD Add random values with a flat distribution to the components of the reciprocal lattice vectors written to the stream, simulating indexing errors. The maximum value that will be added is +/- \fIval\fR percent. -.PD 0 -.B -.IP "\fB--pgraph=\fR\fIfilename\fR -.PD -Write reflection statistics to \fifilename\fR. The file consists of a one-line -header and has columns as follows: the centre of the resolution bin (in inverse -nanometres), the number of reflections in that bin across the whole simulated -data set, the mean partiality for the bin, and the maximum partiality which was -encountered in the bin. - - .SH AUTHOR This page was written by Thomas White. diff --git a/libcrystfel/src/cell.c b/libcrystfel/src/cell.c index 8412fc62..2e87b770 100644 --- a/libcrystfel/src/cell.c +++ b/libcrystfel/src/cell.c @@ -690,7 +690,7 @@ static int same_vector(struct cvec a, struct cvec b) /* Attempt to make 'cell' fit into 'template' somehow */ UnitCell *match_cell(UnitCell *cell, UnitCell *template, int verbose, - float *tols, int reduce) + const float *tols, int reduce) { signed int n1l, n2l, n3l; double asx, asy, asz; diff --git a/libcrystfel/src/cell.h b/libcrystfel/src/cell.h index cbe84772..bd2719dd 100644 --- a/libcrystfel/src/cell.h +++ b/libcrystfel/src/cell.h @@ -118,7 +118,7 @@ extern UnitCell *rotate_cell(UnitCell *in, double omega, double phi, extern void cell_print(UnitCell *cell); extern UnitCell *match_cell(UnitCell *cell, UnitCell *tempcell, int verbose, - float *ltl, int reduce); + const float *ltl, int reduce); extern UnitCell *match_cell_ab(UnitCell *cell, UnitCell *tempcell); diff --git a/libcrystfel/src/index.c b/libcrystfel/src/index.c index 3d0f164b..2baa2644 100644 --- a/libcrystfel/src/index.c +++ b/libcrystfel/src/index.c @@ -157,7 +157,7 @@ void map_all_peaks(struct image *image) void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, int cellr, int verbose, IndexingPrivate **ipriv, - int config_insane, float *ltl) + int config_insane, const float *ltl) { int i; int n = 0; diff --git a/libcrystfel/src/index.h b/libcrystfel/src/index.h index c4e15f09..9d23f3fb 100644 --- a/libcrystfel/src/index.h +++ b/libcrystfel/src/index.h @@ -74,7 +74,8 @@ extern void map_all_peaks(struct image *image); extern void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm, int cellr, int verbose, - IndexingPrivate **priv, int config_insane, float *ltl); + IndexingPrivate **priv, int config_insane, + const float *ltl); extern void cleanup_indexing(IndexingPrivate **priv); diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c index 3e29a4f8..0032a095 100644 --- a/libcrystfel/src/mosflm.c +++ b/libcrystfel/src/mosflm.c @@ -379,7 +379,6 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) rval = read(mosflm->pty, mosflm->rbuffer+mosflm->rbufpos, mosflm->rbuflen-mosflm->rbufpos); - if ( (rval == -1) || (rval == 0) ) return 1; mosflm->rbufpos += rval; @@ -436,7 +435,6 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm) break; case MOSFLM_INPUT_PROMPT : - mosflm_send_next(image, mosflm); endbit_length = i+7; break; @@ -513,6 +511,7 @@ void run_mosflm(struct image *image, UnitCell *cell) remove(mosflm->newmatfile); mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL); + if ( mosflm->pid == -1 ) { ERROR("Failed to fork for MOSFLM\n"); free(mosflm); @@ -545,6 +544,7 @@ void run_mosflm(struct image *image, UnitCell *cell) mosflm->step = 1; /* This starts the "initialisation" procedure */ mosflm->finished_ok = 0; + do { fd_set fds; @@ -559,7 +559,7 @@ void run_mosflm(struct image *image, UnitCell *cell) sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv); - if ( sval == -1 ) { + if ( sval == -1 ) { int err = errno; ERROR("select() failed: %s\n", strerror(err)); rval = 1; @@ -570,6 +570,7 @@ void run_mosflm(struct image *image, UnitCell *cell) rval = 1; } + } while ( !rval ); close(mosflm->pty); diff --git a/libcrystfel/src/peaks.c b/libcrystfel/src/peaks.c index 477d9345..583dff4c 100644 --- a/libcrystfel/src/peaks.c +++ b/libcrystfel/src/peaks.c @@ -148,10 +148,9 @@ static int cull_peaks(struct image *image) /* Returns non-zero if peak has been vetoed. * i.e. don't use result if return value is not zero. */ -static int integrate_peak(struct image *image, int cfs, int css, - double *pfs, double *pss, - double *intensity, double *sigma, - double ir_inn, double ir_mid, double ir_out) +int integrate_peak(struct image *image, int cfs, int css, + double *pfs, double *pss, double *intensity, double *sigma, + double ir_inn, double ir_mid, double ir_out) { signed int fs, ss; double lim_sq, out_lim_sq, mid_lim_sq; diff --git a/src/indexamajig.c b/src/indexamajig.c index 8c073e2c..223403ad 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. * @@ -44,6 +45,8 @@ #include <hdf5.h> #include <gsl/gsl_errno.h> #include <pthread.h> +#include <sys/wait.h> +#include <fcntl.h> #ifdef HAVE_CLOCK_GETTIME #include <time.h> @@ -67,7 +70,6 @@ /* Write statistics at APPROXIMATELY this interval */ #define STATS_EVERY_N_SECONDS (5) - enum { PEAK_ZAEF, PEAK_HDF5, @@ -75,7 +77,7 @@ enum { /* Information about the indexing process which is common to all patterns */ -struct static_index_args +struct index_args { UnitCell *cell; int config_cmfilter; @@ -104,42 +106,23 @@ struct static_index_args double ir_out; /* Output stream */ - pthread_mutex_t *output_mutex; /* Protects the output stream */ FILE *ofh; const struct copy_hdf5_field *copyme; + char *outfile; }; /* Information about the indexing process for one pattern */ -struct index_args +struct pattern_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; -}; - - static void show_help(const char *s) { printf("Syntax: %s [options]\n\n", s); @@ -234,80 +217,108 @@ 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) +static char *get_pattern(FILE *fh, char **use_this_one_instead, + int config_basename, const char *prefix) +{ + char *line; + char *filename; + + /* Get the next filename */ + if ( *use_this_one_instead != NULL ) { + + line = *use_this_one_instead; + *use_this_one_instead = NULL; + + } else { + + char *rval; + + line = malloc(1024*sizeof(char)); + rval = fgets(line, 1023, fh); + if ( rval == NULL ) { + free(line); + return NULL; + } + + } + + chomp(line); + + if ( config_basename ) { + char *tmp; + tmp = safe_basename(line); + free(line); + line = tmp; + } + + filename = malloc(strlen(prefix)+strlen(line)+1); + + snprintf(filename, 1023, "%s%s", prefix, line); + + free(line); + + return filename; +} + + +static void process_image(const struct index_args *iargs, + struct pattern_args *pargs, 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; + 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 r, check; + struct hdfile *hdfile; + char *outfile = iargs->outfile; + struct image image; + char *outfilename = iargs->outfile; + int fd; + FILE *fh; + struct flock fl = {F_WRLCK, SEEK_SET, 0, 0, 0}; image.features = NULL; image.data = NULL; image.flags = NULL; image.indexed_cell = NULL; - image.id = cookie; - image.filename = filename; - image.det = copy_geom(pargs->static_args.det); - image.copyme = pargs->static_args.copyme; + image.det = copy_geom(iargs->det); + image.copyme = iargs->copyme; image.beam = beam; + image.id = cookie; + image.filename = pargs->filename; - pargs->indexable = 0; - - hdfile = hdfile_open(filename); + hdfile = hdfile_open(image.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; - } - + 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); + check = hdf5_read(hdfile, &image, 1); + if ( check ) { + hdfile_close(hdfile); + return; + } - if ( (image.width != image.det->max_fs+1) - || (image.height != image.det->max_ss+1) ) + 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"); + " - 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); + image.width, image.height, + image.det->max_fs + 1, image.det->max_ss + 1); hdfile_close(hdfile); free_detector_geometry(image.det); return; @@ -316,12 +327,12 @@ static void process_image(void *pp, int cookie) if ( image.lambda < 0.0 ) { if ( beam != NULL ) { ERROR("Using nominal photon energy of %.2f eV\n", - beam->photon_energy); + beam->photon_energy); image.lambda = ph_en_to_lambda( - eV_to_J(beam->photon_energy)); + eV_to_J(beam->photon_energy)); } else { ERROR("No wavelength in file, so you need to give " - "a beam parameters file with -b.\n"); + "a beam parameters file with -b.\n"); hdfile_close(hdfile); free_detector_geometry(image.det); return; @@ -335,7 +346,7 @@ static void process_image(void *pp, int cookie) /* Take snapshot of image after CM subtraction but before * the aggressive noise filter. */ - data_size = image.width*image.height*sizeof(float); + data_size = image.width * image.height * sizeof(float); data_for_measurement = malloc(data_size); if ( config_noisefilter ) { @@ -344,25 +355,25 @@ static void process_image(void *pp, int cookie) 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) ) - { + switch ( iargs->peaks ) { + + case PEAK_HDF5: + // Get peaks from HDF5 + if (get_peaks(&image, hdfile, + iargs->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); + case PEAK_ZAEF: + search_peaks(&image, iargs->threshold, + iargs->min_gradient, + iargs->min_snr, + iargs->ir_inn, + iargs->ir_mid, + iargs->ir_out); break; + } /* Get rid of noise-filtered version at this point @@ -374,40 +385,55 @@ static void process_image(void *pp, int cookie) 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 ) { + 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 ) { - + 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); - + iargs->config_closer, + iargs->config_bgsub, + iargs->min_int_snr, + iargs->ir_inn, + iargs->ir_mid, + iargs->ir_out); } - } else { - image.reflections = NULL; + } + /* Write Lock */ + fl.l_pid = getpid(); + + fd = open(outfile, O_WRONLY); + if ( fd == -1) { + ERROR("Error on opening\n"); + exit(1); + } + if (fcntl(fd, F_SETLKW, &fl) == -1) { + ERROR("Error on setting lock wait\n"); + exit(1); } - 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); + fh = fopen(outfilename, "a"); + if ( fh == NULL ) { + ERROR("Couldn't open stream '%s'.\n", outfilename); + } + write_chunk(fh, &image, hdfile, iargs->stream_flags); + fclose(fh); + + /* Unlock stream for other processes */ + fl.l_type = F_UNLCK; /* set to unlock same region */ + if ( fcntl(fd, F_SETLK, &fl) == -1 ) { + ERROR("fcntl"); + exit(1); + } + close(fd); /* Only free cell if found */ cell_free(image.indexed_cell); @@ -421,51 +447,60 @@ static void process_image(void *pp, int cookie) } -static void *get_image(void *qp) +static void run_work(const struct index_args *iargs, + int filename_pipe, int results_pipe, int cookie) { - 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)); + int allDone = 0; + FILE *fh; - /* Get the next filename */ - if ( qargs->use_this_one_instead != NULL ) { + fh = fdopen(filename_pipe, "r"); + if ( fh == NULL ) { + ERROR("Failed to fdopen() the filename pipe!\n"); + close(filename_pipe); + close(results_pipe); + return; + } - line = qargs->use_this_one_instead; - qargs->use_this_one_instead = NULL; + while ( !allDone ) { - } else { + struct pattern_args pargs; + int w, c; + char buf[1024]; + char *line; + char *rval; line = malloc(1024*sizeof(char)); - rval = fgets(line, 1023, qargs->fh); + rval = fgets(line, 1023, fh); if ( rval == NULL ) { - free(pargs); free(line); - return NULL; + if ( feof(fh) ) { + allDone = 1; + continue; + } else { + ERROR("Read error!\n"); + return; + } } - chomp(line); - } - - if ( qargs->config_basename ) { - char *tmp; - tmp = safe_basename(line); - free(line); - line = tmp; - } + chomp(line); + pargs.filename = line; + pargs.indexable = 0; - pargs->filename = malloc(strlen(qargs->prefix)+strlen(line)+1); + process_image(iargs, &pargs, cookie); - snprintf(pargs->filename, 1023, "%s%s", qargs->prefix, line); + /* Request another image */ + c = sprintf(buf, "%i\n", pargs.indexable); + w = write(results_pipe, buf, c); + if ( w < 0 ) { + ERROR("write P0\n"); + } - free(line); + free(line); - return pargs; + } + /* close my pipes */ + fclose(fh); + close(results_pipe); } @@ -492,34 +527,6 @@ static time_t get_monotonic_seconds() #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) @@ -553,7 +560,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; @@ -584,25 +590,31 @@ 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; float ir_inn = 4.0; float ir_mid = 5.0; float ir_out = 7.0; + int n_indexable, n_processed, n_indexable_last_stats; + int n_processed_last_stats; + int t_last_stats; + pid_t *pids; + int *filename_pipes; + FILE **result_fhs; + fd_set fds; + int i; + int allDone; + int *finished; copyme = new_copy_hdf5_field_list(); if ( copyme == NULL ) { @@ -653,7 +665,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:wp:j:x:g:t:o:b:e:", + while ((c = getopt_long(argc, argv, "hi:o:z:p:x:j:g:t:b:e:", longopts, NULL)) != -1) { switch (c) { @@ -683,7 +695,7 @@ int main(int argc, char *argv[]) break; case 'j' : - nthreads = atoi(optarg); + n_proc = atoi(optarg); break; case 'g' : @@ -725,39 +737,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 : @@ -794,12 +777,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("-"); } @@ -826,7 +803,6 @@ int main(int argc, char *argv[]) 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"); @@ -859,8 +835,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; } @@ -966,22 +942,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; @@ -990,7 +964,6 @@ 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; /* Prepare the indexer */ if ( indm != NULL ) { @@ -1006,50 +979,258 @@ 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); + /* 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.ofh = ofh; + 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; + iargs.outfile = outfile; + + n_indexable = 0; + n_processed = 0; + n_indexable_last_stats = 0; + n_processed_last_stats = 0; + t_last_stats = get_monotonic_seconds(); + + FD_ZERO(&fds); + filename_pipes = calloc(n_proc, sizeof(int)); + result_fhs = calloc(n_proc, sizeof(FILE *)); + if ( filename_pipes == NULL ) { + ERROR("Couldn't allocate memory for pipes.\n"); + return 1; + } + if ( result_fhs == NULL ) { + ERROR("Couldn't allocate memory for pipe file handles.\n"); + return 1; + } + + pids = calloc(n_proc, sizeof(pid_t)); + if ( pids == NULL ) { + ERROR("Couldn't allocate memory for PIDs.\n"); + return 1; + } + + finished = calloc(n_proc, sizeof(int)); + if ( finished == NULL ) { + ERROR("Couldn't allocate memory for process flags.\n"); + return 1; + } + + /* Fork the right number of times */ + for ( i=0; i<n_proc; i++ ) { + + pid_t p; + int filename_pipe[2]; + int result_pipe[2]; + + if ( pipe(filename_pipe) == - 1 ) { + ERROR("pipe() failed!\n"); + return 1; + } + + if ( pipe(result_pipe) == - 1 ) { + ERROR("pipe() failed!\n"); + return 1; + } + + p = fork(); + if ( p == -1 ) { + ERROR("fork() failed!\n"); + return 1; + } + + if ( p == 0 ) { + /* Child process gets the 'read' end of the filename + * pipe, and the 'write' end of the result pipe. */ + close(filename_pipe[1]); + close(result_pipe[0]); + run_work(&iargs, filename_pipe[0], result_pipe[1], i); + exit(0); + } + + /* Parent process gets the 'write' end of the filename pipe + * and the 'read' end of the result pipe. */ + pids[i] = p; + close(filename_pipe[0]); + close(result_pipe[1]); + filename_pipes[i] = filename_pipe[1]; + + result_fhs[i] = fdopen(result_pipe[0], "r"); + if ( result_fhs[i] == NULL ) { + ERROR("fdopen() failed.\n"); + return 1; + } + + } + + /* Send first image to all children */ + for ( i=0; i<n_proc; i++ ) { + + char *nextImage; + + nextImage = get_pattern(fh, &use_this_one_instead, + config_basename, prefix); + + if ( nextImage != NULL ) { + + write(filename_pipes[i], nextImage, strlen(nextImage)); + write(filename_pipes[i], "\n", 1); + + free(nextImage); + + } else { + + /* No more files to process.. already? */ + close(filename_pipes[i]); + + } + + } + + allDone = 0; + while ( !allDone ) { + + int r, i; + struct timeval tv; + fd_set fds; + double tNow; + int fdmax; + + tv.tv_sec = 5; + tv.tv_usec = 0; + + FD_ZERO(&fds); + fdmax = 0; + for ( i=0; i<n_proc; i++ ) { + + int fd; + + if ( finished[i] ) continue; + + fd = fileno(result_fhs[i]); + FD_SET(fd, &fds); + if ( fd > fdmax ) fdmax = fd; + + } + + r = select(fdmax+1, &fds, NULL, NULL, &tv); + + if ( r == -1 ) { + ERROR("select() failed: %s\n", strerror(errno)); + continue; + } + + if ( r == 0 ) { + STATUS("Timeout\n"); + continue; + } + + for ( i=0; i<n_proc; i++ ) { + + char *nextImage; + char results[1024]; + char *rval; + int fd; + + if ( finished[i] ) continue; + + fd = fileno(result_fhs[i]); + if ( !FD_ISSET(fd, &fds) ) continue; + + rval = fgets(results, 1024, result_fhs[i]); + if ( rval == NULL ) { + if ( feof(result_fhs[i]) ) { + /* Process died */ + finished[i] = 1; + } else { + ERROR("fgets() failed: %s\n", + strerror(errno)); + } + continue; + } + + chomp(results); + n_indexable += atoi(results); + n_processed++; + + /* Send next filename */ + nextImage = get_pattern(fh, + &use_this_one_instead, + config_basename, + prefix); + + if ( nextImage == NULL ) { + /* No more images */ + finished[i] = 1; + } else { + + r = write(filename_pipes[i], nextImage, + strlen(nextImage)); + r -= write(filename_pipes[i], "\n", 1); + if ( r < 0 ) { + ERROR("write pipe\n"); + } + } + + } + + /* Update progress */ + tNow = get_monotonic_seconds(); + if ( tNow >= t_last_stats+STATS_EVERY_N_SECONDS ) { + + STATUS("%i out of %i indexed so far," + " %i out of %i since the last message.\n", + n_indexable, n_processed, + n_indexable - n_indexable_last_stats, + n_processed - n_processed_last_stats); + + n_indexable_last_stats = n_indexable; + n_processed_last_stats = n_processed; + t_last_stats = tNow; + + } + + allDone = 1; + for ( i=0; i<n_proc; i++ ) { + if ( !finished[i] ) allDone = 0; + } + + } + + for ( i=0; i<n_proc; i++ ) { + close(filename_pipes[i]); + fclose(result_fhs[i]); + } + free(filename_pipes); + free(result_fhs); + free(pids); + free(finished); cleanup_indexing(ipriv); @@ -1062,11 +1243,12 @@ int main(int argc, char *argv[]) free(hdf5_peak_path); free_copy_hdf5_field_list(copyme); cell_free(cell); + free(use_this_one_instead); 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); + n_processed, n_indexable); return 0; } diff --git a/src/partial_sim.c b/src/partial_sim.c index 39621a6c..31822eca 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -178,8 +178,6 @@ static void show_help(const char *s) " -c, --cnoise=<val> Add random noise, with a flat distribution, to the\n" " reciprocal lattice vector components given in the\n" " stream, with maximum error +/- <val> percent.\n" -" --pgraph=<file> Write reflection counts and partiality statistics\n" -" to <file>.\n" "\n" ); } @@ -540,8 +538,6 @@ int main(int argc, char *argv[]) if ( fh != NULL ) { - fprintf(fh, "1/d_nm #refl pmean pmax\n"); - for ( i=0; i<NBINS; i++ ) { double rcen; diff --git a/tests/symmetry_check.c b/tests/symmetry_check.c index 7c06dfc8..edbacbf1 100644 --- a/tests/symmetry_check.c +++ b/tests/symmetry_check.c @@ -324,11 +324,6 @@ int main(int argc, char *argv[]) check_subgroup("4/m", "-4", 1, 1, 2, &fail); check_subgroup("622", "321_H", 1, 1, 2, &fail); - check_subgroup("4/mmm", "-42m", 1, 1, 2, &fail); - check_subgroup("4/mmm", "-4m2", 1, 1, 2, &fail); - check_subgroup("4/mmm", "4/m", 1, 1, 2, &fail); - check_subgroup("4/mmm", "4mm", 1, 1, 2, &fail); - /* Tetartohedral */ check_subgroup("6/mmm", "-3_H", 1, 1, 4, &fail); |