From 7ce60a63de7750123d94ac77318547ba6e8b4cf1 Mon Sep 17 00:00:00 2001 From: Chuck Date: Thu, 29 Mar 2012 12:22:06 +0200 Subject: Fork indexamajig --- src/indexamajig.c | 880 +++++++++++++++++++++++++++++++++--------------------- 1 file changed, 534 insertions(+), 346 deletions(-) (limited to 'src') diff --git a/src/indexamajig.c b/src/indexamajig.c index efbf4b44..9ddbb6ea 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -1,32 +1,12 @@ /* - * indexamajig.c + * indexamajigFork.c * * Index patterns, output hkl+intensity etc. * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. - * Copyright © 2012 Richard Kirian - * Copyright © 2012 Lorenzo Galli + * (c) 2006-2012 Thomas White + * (c) 2012- Chunhong Yoon * - * Authors: - * 2010-2012 Thomas White - * 2011 Richard Kirian - * 2012 Lorenzo Galli - * - * This file is part of CrystFEL. - * - * CrystFEL is free software: you can redistribute it and/or modify - * it under the terms of the GNU General Public License as published by - * the Free Software Foundation, either version 3 of the License, or - * (at your option) any later version. - * - * CrystFEL is distributed in the hope that it will be useful, - * but WITHOUT ANY WARRANTY; without even the implied warranty of - * MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the - * GNU General Public License for more details. - * - * You should have received a copy of the GNU General Public License - * along with CrystFEL. If not, see . + * Part of CrystFEL - crystallography with a FEL * */ @@ -44,6 +24,8 @@ #include #include #include +#include +#include #ifdef HAVE_CLOCK_GETTIME #include @@ -51,22 +33,46 @@ #include #endif -#include "utils.h" -#include "hdf5-file.h" -#include "index.h" -#include "peaks.h" -#include "detector.h" -#include "filters.h" -#include "thread-pool.h" -#include "beam-parameters.h" -#include "geometry.h" -#include "stream.h" -#include "reflist-utils.h" +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include +#include /* Write statistics at APPROXIMATELY this interval */ -#define STATS_EVERY_N_SECONDS (5) +#define STATS_EVERY_N_SECONDS (2) + +#define LINE_LENGTH 1024 +#define BUFSIZE 64 + +#ifdef HAVE_CLOCK_GETTIME +static double get_time() +{ + struct timespec tp; + clock_gettime(CLOCK_MONOTONIC, &tp); + double sec = (double) tp.tv_sec+ (double) tp.tv_nsec/1000000000; + return sec; //nano resolution +} +#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 double get_time() +{ + struct timeval tp; + gettimeofday(&tp, NULL); + double sec = (double) tp.tv_sec+ (double) tp.tv_usec/1000000; + return sec; //micro resolution +} +#endif enum { PEAK_ZAEF, @@ -82,6 +88,7 @@ struct static_index_args int config_noisefilter; int config_verbose; int stream_flags; /* What goes into the output? */ + int config_polar; int config_satcorr; int config_closer; int config_insane; @@ -102,11 +109,14 @@ struct static_index_args double ir_inn; double ir_mid; double ir_out; - + /* Output stream */ pthread_mutex_t *output_mutex; /* Protects the output stream */ FILE *ofh; + + char *outfile; const struct copy_hdf5_field *copyme; + int nProcesses; }; @@ -136,9 +146,45 @@ struct queue_args int n_processed; int n_indexable_last_stats; int n_processed_last_stats; - int t_last_stats; + + int n_indexableTotal; + int n_processedTotal; + int n_indexable_last_statsTotal; + int n_processed_last_statsTotal; + + int nPerProcess; + int updateReader; + + double t_last_stats; + }; +// Count number of patterns in .lst +int count_patterns(FILE *fh){ + char *rval; + int n_total = 0; + do{ + char line[LINE_LENGTH]; + rval = fgets(line,LINE_LENGTH-1,fh); + if (rval != NULL){ + n_total++; + } + }while(rval!=NULL); + if (ferror(fh)) { + printf("Read error\n"); + return -1; + } + return n_total; +} + +// Assign a batch number to a process +int getBatchNum(int pid[], int length){ + int i, id = 0; + for (i=0; i0)*pow(2,i)); // child = 1, parent = 0 + } + return id; +} static void show_help(const char *s) { @@ -189,34 +235,35 @@ static void show_help(const char *s) "The default is '--record=integrated'.\n" "\n\n" "For more control over the process, you might need:\n\n" -" --cell-reduction= Use as the cell reduction method. Choose from:\n" -" none : no matching, just use the raw cell.\n" -" reduce : full cell reduction.\n" -" compare : match by at most changing the order of\n" -" the indices.\n" -" compare_ab : compare 'a' and 'b' lengths only.\n" -" --tolerance= Set the tolerances for cell reduction.\n" -" Default: 5,5,5,1.5.\n" -" --filter-cm Perform common-mode noise subtraction on images\n" -" before proceeding. Intensities will be extracted\n" -" from the image as it is after this processing.\n" -" --filter-noise Apply an aggressive noise filter which sets all\n" -" pixels in each 3x3 region to zero if any of them\n" -" have negative values. Intensity measurement will\n" -" be performed on the image as it was before this.\n" -" --no-sat-corr Don't correct values of saturated peaks using a\n" -" table included in the HDF5 file.\n" -" --threshold= Only accept peaks above ADU. Default: 800.\n" -" --min-gradient= Minimum gradient for Zaefferer peak search.\n" -" Default: 100,000.\n" -" --min-snr= Minimum signal-to-noise ratio for peaks.\n" -" Default: 5.\n" -" --min-integration-snr= Minimum signal-to-noise ratio for peaks\n" -" during integration. Default: -infinity.\n" -" --int-radius= Set the integration radii. Default: 4,5,7.\n" -"-e, --image= Use this image from the HDF5 file.\n" -" Example: /data/data0.\n" -" Default: The first one found.\n" +" --cell-reduction= Use as the cell reduction method. Choose from:\n" +" none : no matching, just use the raw cell.\n" +" reduce : full cell reduction.\n" +" compare : match by at most changing the order of\n" +" the indices.\n" +" compare_ab : compare 'a' and 'b' lengths only.\n" +" --tolerance= Set the tolerance for a,b,c axis (in %%)\n" +" and for the angles (in deg) when reducing\n" +" or comparing (default is 5%% and 1.5deg)\n" +" --filter-cm Perform common-mode noise subtraction on images\n" +" before proceeding. Intensities will be extracted\n" +" from the image as it is after this processing.\n" +" --filter-noise Apply an aggressive noise filter which sets all\n" +" pixels in each 3x3 region to zero if any of them\n" +" have negative values. Intensity measurement will\n" +" be performed on the image as it was before this.\n" +" --unpolarized Don't correct for the polarisation of the X-rays.\n" +" --no-sat-corr Don't correct values of saturated peaks using a\n" +" table included in the HDF5 file.\n" +" --threshold= Only accept peaks above ADU. Default: 800.\n" +" --min-gradient= Minimum gradient for Zaefferer peak search.\n" +" Default: 100,000.\n" +" --min-snr= Minimum signal-to-noise ratio for peaks.\n" +" Default: 5.\n" +" --min-integration-snr= Minimum signal-to-noise ratio for peaks\n" +" during integration. Default: -infinity.\n" +" -e, --image= Use this image from the HDF5 file.\n" +" Example: /data/data0.\n" +" Default: The first one found.\n" "\n" "\nFor time-resolved stuff, you might want to use:\n\n" " --copy-hdf5-field Copy the value of field into the stream. You\n" @@ -245,289 +292,291 @@ static void show_help(const char *s) ); } +int readUpdate(void *qp, int fd_pipe[][2], int finish){ + struct queue_args *qargs = qp; + int i,j; + int rFlag; + int pipesOpen = 0; + char bufferR[BUFSIZE]; + memset(bufferR, 0, BUFSIZE); + int numFields = 2; + + for (i=0; istatic_args.nProcesses; i++){ + if (finish && i == qargs->updateReader) { + // do nothing + }else{ + rFlag = read(fd_pipe[i][0], bufferR, BUFSIZE-1); // wait till something to read + if (rFlag != 0){ + char * pch; + char delims[] = "#"; + pch = strtok (bufferR, delims); + for (j=0;jn_indexableTotal += atoi(pch); + break; + case 1: + qargs->n_processedTotal += atoi(pch); + break; + } + pch = strtok (NULL, delims); + } + memset(bufferR, 0, BUFSIZE); + pipesOpen++; + } + } + } + STATUS("%i out of %i indexed so far," + " %i out of %i since the last message.\n\n", + qargs->n_indexableTotal, qargs->n_processedTotal, + qargs->n_indexableTotal-qargs->n_indexable_last_statsTotal, + qargs->n_processedTotal-qargs->n_processed_last_statsTotal); -static void process_image(void *pp, int cookie) -{ - struct index_args *pargs = pp; - struct hdfile *hdfile; - struct image image; + qargs->n_indexable_last_statsTotal = qargs->n_indexableTotal; + qargs->n_processed_last_statsTotal = qargs->n_processedTotal; + + return pipesOpen; +} + +// Manipulate image +static void process_image(char ***array, int batch, void *qp, int fd_pipe[][2]) { + + struct queue_args *qargs = qp; 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.id = cookie; - image.filename = filename; - image.det = copy_geom(pargs->static_args.det); - image.copyme = pargs->static_args.copyme; - image.beam = beam; - - 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"); - } - - pargs->indexable = 0; + UnitCell *cell = qargs->static_args.cell; + int config_cmfilter = qargs->static_args.config_cmfilter; + int config_noisefilter = qargs->static_args.config_noisefilter; + int config_verbose = qargs->static_args.config_verbose; + //int config_polar = qargs->static_args.config_polar; + IndexingMethod *indm = qargs->static_args.indm; + struct beam_params *beam = qargs->static_args.beam; + + int i; + int r, check; + struct hdfile *hdfile; + char *outfile = qargs->static_args.outfile; + int nPerBatch = qargs->nPerProcess; + + for (i=0; istatic_args.det); + image.copyme = qargs->static_args.copyme; + image.beam = beam; + image.id = batch; // MUST SET ID FOR MOSFLM TO WORK PROPERLY + + 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"); + } - hdfile = hdfile_open(filename); - if ( hdfile == NULL ) return; + //pargs->indexable = 0; - if ( pargs->static_args.element != NULL ) { + char *filename = NULL; + char *line = array[batch][i]; + chomp(line); +//printf("%d ***************%s\n",batch,line); + filename = malloc(strlen(qargs->prefix)+strlen(line)+1); + snprintf(filename,LINE_LENGTH-1,"%s%s",qargs->prefix,line); //maximum print length + free(line); + image.filename = filename; +//printf("%d ***************%s\n",batch,filename); + hdfile = hdfile_open(filename); + if (hdfile == NULL) return; - int r; - r = hdfile_set_image(hdfile, pargs->static_args.element); - if ( r ) { - ERROR("Couldn't select path '%s'\n", - pargs->static_args.element); + r = hdfile_set_first_image(hdfile, "/"); // Need this to read hdf5 files + if (r) { + ERROR("Couldn't select first path\n"); hdfile_close(hdfile); return; } - } else { - - int r; - r = hdfile_set_first_image(hdfile, "/"); - if ( r ) { - ERROR("Couldn't select first path\n"); + check = hdf5_read(hdfile,&image,1); + if (check == 1){ 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"); + 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; } - } - fill_in_values(image.det, hdfile); - if ( config_cmfilter ) { - filter_cm(&image); - } + 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); - /* 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_cmfilter ) { + filter_cm(&image); + } - if ( config_noisefilter ) { - filter_noise(&image, data_for_measurement); - } else { - memcpy(data_for_measurement, image.data, data_size); - } + /* 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); - 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"); + if ( config_noisefilter ) { + filter_noise(&image, data_for_measurement); + } else { + memcpy(data_for_measurement, image.data, data_size); } - break; + switch ( qargs->static_args.peaks ) + { + case PEAK_HDF5 : + /* Get peaks from HDF5 */ + if ( get_peaks(&image, hdfile, + qargs->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, + search_peaks(&image, qargs->static_args.threshold, + qargs->static_args.min_gradient, + qargs->static_args.min_snr, + qargs->static_args.ir_inn, + qargs->static_args.ir_mid, + qargs->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; + + ///// RUN INDEXING HERE ///// + index_pattern(&image, cell, indm, qargs->static_args.cellr, + config_verbose, qargs->static_args.ipriv, + qargs->static_args.config_insane, qargs->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); - + if ( image.reflections != NULL ) { + integrate_reflections(&image, + qargs->static_args.config_closer, + qargs->static_args.config_bgsub, + qargs->static_args.min_int_snr, + qargs->static_args.ir_inn, + qargs->static_args.ir_mid, + qargs->static_args.ir_out); + } + } else { + image.reflections = NULL; } - } 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; + // Write Lock + struct flock fl = {F_WRLCK, SEEK_SET, 0, 0, 0 }; + int fd; + fl.l_pid = getpid(); + + char *outfilename = NULL; + chomp(outfile); // may not need this +//printf("%d ***************%s\n",batch,outfile); + outfilename = malloc(strlen(outfile)+1); + snprintf(outfilename,LINE_LENGTH-1,"%s",outfile); //maximum print length +//printf("%d ***************%s\n",batch,outfilename); + if ((fd = open(outfilename, O_WRONLY)) == -1) { + perror("Error on opening\n"); + exit(1); + } + + if (fcntl(fd, F_SETLKW, &fl) == -1) { + perror("Error on setting lock wait\n"); + exit(1); + } + + // LOCKED! Write chunk + FILE *fh; +//printf("%d ***************%s\n",batch,outfilename); + fh = fopen(outfilename,"a"); + if (fh == NULL) { + perror("Error inside lock\n"); + } + write_chunk(fh, &image, hdfile, qargs->static_args.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) { + perror("fcntl"); + exit(1); + } + close(fd); + + ///// WRITE UPDATE ///// + double seconds; + qargs->n_indexable += ( image.indexed_cell != NULL ); + qargs->n_processed++; + seconds = get_time(); + if ( seconds >= qargs->t_last_stats+STATS_EVERY_N_SECONDS + || qargs->n_processed == qargs->nPerProcess) { // Must write if finished + + // WRITE PIPE HERE + char bufferW[BUFSIZE]; + memset(bufferW, 0, BUFSIZE); + sprintf(bufferW,"%d#%d", + qargs->n_indexable - qargs->n_indexable_last_stats, + qargs->n_processed - qargs->n_processed_last_stats); + write(fd_pipe[batch][1], bufferW, BUFSIZE-1); + + // Update stats + qargs->n_processed_last_stats = qargs->n_processed; + qargs->n_indexable_last_stats = qargs->n_indexable; + qargs->t_last_stats = seconds; + + ///// READ UPDATE ///// + if (batch == qargs->updateReader) { + readUpdate(qargs,fd_pipe, 0); + } } - 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 ///// + /* 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); } - - free(pargs->filename); - free(pargs); } - static int parse_cell_reduction(const char *scellr, int *err, int *reduction_needs_cell) { @@ -565,6 +614,7 @@ int main(int argc, char *argv[]) int config_cmfilter = 0; int config_noisefilter = 0; int config_verbose = 0; + int config_polar = 1; int config_satcorr = 1; int config_checkprefix = 1; int config_closer = 1; @@ -591,7 +641,7 @@ 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; + int nProcesses = 1; pthread_mutex_t output_mutex = PTHREAD_MUTEX_INITIALIZER; char *prepare_line; char prepare_filename[1024]; @@ -610,7 +660,7 @@ int main(int argc, char *argv[]) float ir_inn = 4.0; float ir_mid = 5.0; float ir_out = 7.0; - + copyme = new_copy_hdf5_field_list(); if ( copyme == NULL ) { ERROR("Couldn't allocate HDF5 field list.\n"); @@ -689,7 +739,7 @@ int main(int argc, char *argv[]) break; case 'j' : - nthreads = atoi(optarg); + nProcesses = atoi(optarg); break; case 'g' : @@ -828,11 +878,11 @@ int main(int argc, char *argv[]) } else { ofh = fopen(outfile, "w"); } +//printf("***************%s\n",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"); @@ -865,7 +915,7 @@ int main(int argc, char *argv[]) } } - if ( nthreads == 0 ) { + if ( nProcesses == 0 ) { ERROR("Invalid number of threads.\n"); return 1; } @@ -964,7 +1014,7 @@ int main(int argc, char *argv[]) cell = NULL; } free(pdb); - + write_stream_header(ofh, argc, argv); if ( beam != NULL ) { @@ -1010,6 +1060,7 @@ int main(int argc, char *argv[]) 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_polar = config_polar; qargs.static_args.config_satcorr = config_satcorr; qargs.static_args.config_closer = config_closer; qargs.static_args.config_insane = config_insane; @@ -1034,6 +1085,7 @@ int main(int argc, char *argv[]) qargs.static_args.stream_flags = stream_flags; qargs.static_args.hdf5_peak_path = hdf5_peak_path; qargs.static_args.copyme = copyme; + qargs.static_args.nProcesses = nProcesses; qargs.static_args.ir_inn = ir_inn; qargs.static_args.ir_mid = ir_mid; qargs.static_args.ir_out = ir_out; @@ -1045,28 +1097,164 @@ int main(int argc, char *argv[]) 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); - 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; + qargs.n_indexableTotal = 0; + qargs.n_processedTotal = 0; + qargs.n_indexable_last_statsTotal = 0; + qargs.n_processed_last_statsTotal = 0; + qargs.t_last_stats = get_time(); + qargs.updateReader = nProcesses-1; // last process reads updates + + //////////// Read .lst file /////////////// + int i,j; + double t1, t2; + rewind(fh); // make sure count from start + n_images = count_patterns(fh); + rewind(fh); + // Divide images into nProcesses + int nPerBatch = (int)ceil((float)n_images/nProcesses); + // Malloc 3D array + char ***array; + int nChars = LINE_LENGTH; + array = malloc(nProcesses * sizeof(char **)); + if (array == NULL){ + printf("Error\n"); + return -1; + } + for(i=0;i= nProcesses){ + exit(0); // kill + } + batchNum = (int) num; + + // Calculate how many images to process + qargs.nPerProcess = 0; + i=batchNum; + for (j=0; j= qargs.t_last_stats+STATS_EVERY_N_SECONDS ) { + pipesOpen = readUpdate(&qargs,fd_pipe,1); + if (pipesOpen == 0) { + // Close all read pipes + for (i=0;i