aboutsummaryrefslogtreecommitdiff
path: root/src/indexamajig.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/indexamajig.c')
-rw-r--r--src/indexamajig.c880
1 files changed, 534 insertions, 346 deletions
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 <taw@physics.org>
+ * (c) 2012- Chunhong Yoon <chun.hong.yoon@desy.de>
*
- * Authors:
- * 2010-2012 Thomas White <taw@physics.org>
- * 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 <http://www.gnu.org/licenses/>.
+ * Part of CrystFEL - crystallography with a FEL
*
*/
@@ -44,6 +24,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>
@@ -51,22 +33,46 @@
#include <sys/time.h>
#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 <crystfel/utils.h>
+#include <crystfel/hdf5-file.h>
+#include <crystfel/index.h>
+#include <crystfel/peaks.h>
+#include <crystfel/detector.h>
+#include <crystfel/filters.h>
+#include <crystfel/thread-pool.h>
+#include <crystfel/beam-parameters.h>
+#include <crystfel/geometry.h>
+#include <crystfel/stream.h>
+#include <crystfel/reflist-utils.h>
/* 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; i<length; i++) {
+ id += ((pid[i]>0)*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=<m> Use <m> 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=<tol> 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=<n> Only accept peaks above <n> ADU. Default: 800.\n"
-" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n"
-" Default: 100,000.\n"
-" --min-snr=<n> Minimum signal-to-noise ratio for peaks.\n"
-" Default: 5.\n"
-" --min-integration-snr=<n> Minimum signal-to-noise ratio for peaks\n"
-" during integration. Default: -infinity.\n"
-" --int-radius=<r> Set the integration radii. Default: 4,5,7.\n"
-"-e, --image=<element> Use this image from the HDF5 file.\n"
-" Example: /data/data0.\n"
-" Default: The first one found.\n"
+" --cell-reduction=<m> Use <m> 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=<a,b,c,angl> 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=<n> Only accept peaks above <n> ADU. Default: 800.\n"
+" --min-gradient=<n> Minimum gradient for Zaefferer peak search.\n"
+" Default: 100,000.\n"
+" --min-snr=<n> Minimum signal-to-noise ratio for peaks.\n"
+" Default: 5.\n"
+" --min-integration-snr=<n> Minimum signal-to-noise ratio for peaks\n"
+" during integration. Default: -infinity.\n"
+" -e, --image=<element> 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 <f> Copy the value of field <f> 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; i<qargs->static_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;j<numFields;j++)
+ {
+ switch (j){
+ case 0:
+ qargs->n_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; i<nPerBatch; i++) {
+ struct image image;
+ image.features = NULL;
+ image.data = NULL;
+ image.flags = NULL;
+ image.indexed_cell = NULL;
+ image.det = copy_geom(qargs->static_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;i++){
+ array[i] = malloc(nPerBatch * sizeof(char *));
+ if(array[i] == NULL){
+ printf("Error\n");
+ return -1;
+ }
+ for(j=0;j<nPerBatch;j++){
+ array[i][j] = malloc(nChars * sizeof(char));
+ if(array[i][j] == NULL){
+ printf("Error\n");
+ return -1;
+ }
+ }
+ }
+ // Read from file
+ for (i=0; i<nProcesses; i++) {
+ for (j=0; j<nPerBatch; j++) {
+ if(fgets(&array[i][j][0],(nChars * sizeof(char))-1,fh) == NULL){
+ array[i][j] = NULL;
+ }
+ }
+ }
+ fclose(fh);
+ printf("\n"); // flush buffer
+
+ // Clear output file content
+ char *myOutfilename = NULL;
+ chomp(prefix);
+ chomp(outfile);
+ myOutfilename = malloc(strlen(outfile)+1);
+ snprintf(myOutfilename,LINE_LENGTH-1,"%s",outfile); //maximum print length
+ FILE *tfh;
+ tfh = fopen(myOutfilename,"a+");
+ if (tfh == NULL) { printf("Error\n"); }
+ fclose(tfh);
+ qargs.static_args.outfile = outfile;
+
+ //////////// PIPING ///////////////
+ char bufferW[BUFSIZE], bufferR[BUFSIZE]; // write and read buffers
+ int fd_pipe[nProcesses][2];
+
+ memset(bufferW, 0, BUFSIZE);
+ memset(bufferR, 0, BUFSIZE);
+
+ for (i=0;i<nProcesses;i++){ // pipe nthread times
+ pipe(fd_pipe[i]);
+ }
+
+ //////////// FORKING ///////////////
+ int power = 10; // 2^power must be larger than nProcesses
+ int pid[power];
+ double num = 0;
+ int batchNum = 0;
+ // Fork 2^power times
+ for (i=0;i<power;i++){
+ pid[i] = fork();
+ }
+ // Assign id
+ for (i=0;i<power;i++){
+ if (pid[i] == 0){ // keep parents and kill off children
+ num += pow(2,i);
+ }
+ }
+ // Kill if batchNum too high
+ if (num >= nProcesses){
+ exit(0); // kill
+ }
+ batchNum = (int) num;
+
+ // Calculate how many images to process
+ qargs.nPerProcess = 0;
+ i=batchNum;
+ for (j=0; j<nPerBatch; j++) {
+ if (array[i][j] != '\0'){
+ qargs.nPerProcess++;
+ }
+ }
+ //printf("%d: %d images\n",batchNum,qargs.nPerProcess);
+
+ ////////// PLUMBING //////////
+ if (batchNum == qargs.updateReader){
+ for (i=0;i<nProcesses;i++){
+ if (i != batchNum) {
+ close(fd_pipe[i][1]); // only open own write pipe
+ }
+ // read from all processes
+ }
+ } else {
+ for (i=0;i<nProcesses;i++){
+ if (i != batchNum) {
+ close(fd_pipe[i][1]); // only open own write pipe
+ }
+ close(fd_pipe[i][0]); // close all read pipes
+ }
+ }
+
+ //////////// INDEXING ///////////////
+ //printf("PID&PPID: %d %d\n",getpid(),getppid());
+
+ t1 = get_time();
+
+ process_image(array,batchNum,&qargs,fd_pipe);
+
+ ///// READ UPDATE /////
+ if (batchNum == qargs.updateReader){ // loop until every process is finished
+ int pipesOpen;
+ int loopFlag = 1;
+ while (loopFlag) {
+ double seconds;
+ seconds = get_time();
+ if ( seconds >= 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<qargs.static_args.nProcesses;i++){
+ close(fd_pipe[i][0]);
+ }
+ loopFlag = 0; // Must execute this to get out
+ }
+ }
+ }
+ }
+
+ close(fd_pipe[batchNum][1]); // close own write pipe
+
+ t2 = get_time();
+
+ //////////// EXITING ///////////////
+
+ printf("%d: Fork Time %.2fs\n",batchNum,t2-t1);
+
+ if (batchNum == qargs.updateReader){
+ STATUS("There were %i images, of which %i could be indexed.\n",
+ n_images, qargs.n_indexableTotal);
+ }
+
+ exit(0);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-17 06:31:40 +0000