aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--.gitignore1
-rw-r--r--libcrystfel/src/index.c42
-rw-r--r--libcrystfel/src/mosflm.c83
-rw-r--r--src/indexamajig.c550
4 files changed, 460 insertions, 216 deletions
diff --git a/.gitignore b/.gitignore
index b8a4028d..fba1c515 100644
--- a/.gitignore
+++ b/.gitignore
@@ -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/libcrystfel/src/index.c b/libcrystfel/src/index.c
index 3d0f164b..5b46edcb 100644
--- a/libcrystfel/src/index.c
+++ b/libcrystfel/src/index.c
@@ -50,6 +50,35 @@
#include "reax.h"
#include "geometry.h"
+#ifdef HAVE_CLOCK_GETTIME
+#include <time.h>
+#else
+#include <sys/time.h>
+#endif
+#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
+
+
+
+
/* Base class constructor for unspecialised indexing private data */
IndexingPrivate *indexing_private(IndexingMethod indm)
@@ -159,6 +188,8 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm,
int cellr, int verbose, IndexingPrivate **ipriv,
int config_insane, float *ltl)
{
+//double t1,t2;
+//t1 = get_time();
int i;
int n = 0;
@@ -167,6 +198,9 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm,
map_all_peaks(image);
image->indexed_cell = NULL;
+double tDI,tDO;
+double tMI,tMO;
+
while ( indm[n] != INDEXING_NONE ) {
image->ncells = 0;
@@ -176,10 +210,16 @@ void index_pattern(struct image *image, UnitCell *cell, IndexingMethod *indm,
case INDEXING_NONE :
return;
case INDEXING_DIRAX :
+tDI = get_time();
run_dirax(image);
+tDO = get_time();
+ERROR("run_dirax DONE %.2f %.2f %.2f\n",tDI,tDO,tDO-tDI);
break;
case INDEXING_MOSFLM :
+tMI = get_time();
run_mosflm(image, cell);
+tMO = get_time();
+ERROR("run_mosflm DONE %.2f %.2f %.5f\n",tMI,tMO,tMO-tMI);
break;
case INDEXING_REAX :
reax_index(ipriv[n], image, cell);
@@ -251,6 +291,8 @@ done:
/* May free(NULL) if all algorithms were tried and no success */
cell_free(image->candidate_cells[i]);
}
+ //t2 = get_time();
+ //ERROR("index_pattern DONE %.2f %.2f %.5f\n",t1,t2,t2-t1);
}
diff --git a/libcrystfel/src/mosflm.c b/libcrystfel/src/mosflm.c
index 3e29a4f8..03e6cbeb 100644
--- a/libcrystfel/src/mosflm.c
+++ b/libcrystfel/src/mosflm.c
@@ -82,6 +82,37 @@
#include "peaks.h"
+
+
+#ifdef HAVE_CLOCK_GETTIME
+#include <time.h>
+#else
+#include <sys/time.h>
+#endif
+#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
+
+
+
+
#define MOSFLM_VERBOSE 0
@@ -358,6 +389,7 @@ static void mosflm_send_next(struct image *image, struct mosflm_data *mosflm)
break;
case 10 :
+//ERROR(">>>>>>>>>>>>>> OK\n");
mosflm_sendline("GO\n", mosflm);
mosflm->finished_ok = 1;
break;
@@ -376,12 +408,12 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm)
{
int rval;
int no_string = 0;
-
+//ERROR("read pty: %d %d %d ===",mosflm->pty,mosflm->rbuffer+mosflm->rbufpos, (int) mosflm->rbuflen-mosflm->rbufpos);
rval = read(mosflm->pty, mosflm->rbuffer+mosflm->rbufpos,
mosflm->rbuflen-mosflm->rbufpos);
-
+//ERROR("read val: %d\n",rval);
if ( (rval == -1) || (rval == 0) ) return 1;
-
+//ERROR("Enter here\n");
mosflm->rbufpos += rval;
assert(mosflm->rbufpos <= mosflm->rbuflen);
@@ -436,7 +468,7 @@ static int mosflm_readable(struct image *image, struct mosflm_data *mosflm)
break;
case MOSFLM_INPUT_PROMPT :
-
+//ERROR("Enter mosflm_send_next\n");
mosflm_send_next(image, mosflm);
endbit_length = i+7;
break;
@@ -512,14 +544,16 @@ void run_mosflm(struct image *image, UnitCell *cell)
snprintf(mosflm->newmatfile, 127, "xfel-%i.newmat", image->id);
remove(mosflm->newmatfile);
+ // fork a new process operating in pseudoterminal
mosflm->pid = forkpty(&mosflm->pty, NULL, NULL, NULL);
+//ERROR("forkpty: %d\n",mosflm->pid);
if ( mosflm->pid == -1 ) {
ERROR("Failed to fork for MOSFLM\n");
free(mosflm);
return;
}
- if ( mosflm->pid == 0 ) {
-
+ if ( mosflm->pid == 0 ) { // child process
+//ERROR("Enter the dragon\n");
/* Child process: invoke MOSFLM */
struct termios t;
@@ -545,32 +579,61 @@ void run_mosflm(struct image *image, UnitCell *cell)
mosflm->step = 1; /* This starts the "initialisation" procedure */
mosflm->finished_ok = 0;
+
+//double t1,t2,t3,t4,t5,t6,t7;
+
+//t1 = get_time();
+
do {
+//t2 = get_time();
+
fd_set fds;
struct timeval tv;
int sval;
+//t3 = get_time();
+
FD_ZERO(&fds);
- FD_SET(mosflm->pty, &fds);
+ FD_SET(mosflm->pty, &fds); // file descriptor set
+
+//t4 = get_time();
- tv.tv_sec = 30;
+ tv.tv_sec = 30; // 30 second timeout
tv.tv_usec = 0;
- sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv);
+//t5 = get_time();
+ //sval = 1;
+ sval = select(mosflm->pty+1, &fds, NULL, NULL, &tv); // is mosflm ready for reading?
+
+ if (sval != 1){
+ ERROR("******************* sval: %d ",sval);
+ }
+//t6 = get_time();
if ( sval == -1 ) {
int err = errno;
ERROR("select() failed: %s\n", strerror(err));
rval = 1;
} else if ( sval != 0 ) {
- rval = mosflm_readable(image, mosflm);
+//ERROR("mosflm_readable: %d %d %d ===",mosflm->pty,mosflm->rbuffer+mosflm->rbufpos, (int) mosflm->rbuflen-mosflm->rbufpos);
+ rval = mosflm_readable(image, mosflm); // read mosflm results
} else {
ERROR("No response from MOSFLM..\n");
rval = 1;
}
+ //if (sval != 1){
+ // ERROR("rval: %d\n",rval);
+ //}
+
+//t7 = get_time();
+//ERROR("mosflm_do DONE %.5f %.5f %.5f %.5f %.5f %.5f\n",t2-t1,t3-t2,t4-t3,t5-t4,t6-t5,t7-t6);
} while ( !rval );
+ //ERROR("DONE\n");
+
+
+
close(mosflm->pty);
free(mosflm->rbuffer);
diff --git a/src/indexamajig.c b/src/indexamajig.c
index 281cac43..9372c447 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>
@@ -51,22 +54,25 @@
#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 LINE_LENGTH 1024
+
+#define BUFFER PIPE_BUF
enum {
PEAK_ZAEF,
@@ -107,6 +113,7 @@ struct static_index_args
pthread_mutex_t *output_mutex; /* Protects the output stream */
FILE *ofh;
const struct copy_hdf5_field *copyme;
+ char *outfile;
};
@@ -137,6 +144,7 @@ struct queue_args
int n_indexable_last_stats;
int n_processed_last_stats;
int t_last_stats;
+ int updateReader;
};
@@ -228,9 +236,8 @@ static void show_help(const char *s)
"\n"
"\nOptions you probably won't need:\n\n"
" --no-check-prefix Don't attempt to correct the --prefix.\n"
-" --no-closer-peak Don't integrate from the location of a nearby peak\n"
-" instead of the position closest to the reciprocal\n"
-" lattice point.\n"
+" --closer-peak Don't integrate from the location of a nearby peak\n"
+" instead of the predicted spot. Don't use.\n"
" --insane Don't check that the reduced cell accounts for at\n"
" least 10%% of the located peaks.\n"
" --no-bg-sub Don't subtract local background estimates from\n"
@@ -246,83 +253,94 @@ static void show_help(const char *s)
}
-static void process_image(void *pp, int cookie)
+// Get next pattern in .lst
+char* get_pattern(FILE *fh) {
+ char *rval;
+ char line[LINE_LENGTH];
+ rval = fgets(line, LINE_LENGTH - 1, fh);
+ if (ferror(fh)) {
+ printf("Read error\n");
+ rval = NULL;
+ }
+ return rval;
+}
+
+
+static void process_image(void *qp, void *pp, int cookie)
{
struct index_args *pargs = pp;
- struct hdfile *hdfile;
- struct image image;
+ 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;
+ 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;
+ IndexingMethod *indm = qargs->static_args.indm;
+ struct beam_params *beam = qargs->static_args.beam;
+ int r, check;
+ struct hdfile *hdfile;
+ char *outfile = qargs->static_args.outfile;
+ struct image image;
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(qargs->static_args.det);
+ image.copyme = qargs->static_args.copyme;
image.beam = beam;
+ image.id = cookie; // MUST SET ID FOR MOSFLM TO WORK PROPERLY
- pargs->indexable = 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");
+ }
+ char *filename = NULL;
+ char *imagename = pargs->filename;
+ chomp(imagename);
+ filename = malloc(strlen(qargs->prefix) + strlen(imagename) + 1);
+ snprintf(filename, LINE_LENGTH - 1, "%s%s", qargs->prefix, imagename);
+ image.filename = filename;
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;
- }
+ if (hdfile == NULL) return;
+ 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;
}
- hdf5_read(hdfile, &image, pargs->static_args.config_satcorr);
+ check = hdf5_read(hdfile, &image, 1);
+ if (check == 1) {
+ 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;
}
- if ( image.lambda < 0.0 ) {
- if ( beam != NULL ) {
+ 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;
@@ -330,40 +348,37 @@ static void process_image(void *pp, int cookie)
}
fill_in_values(image.det, hdfile);
- if ( config_cmfilter ) {
+ 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);
+ // Take snapshot of image after CM subtraction but before
+ // the aggressive noise filter.
+ data_size = image.width * image.height * sizeof (float);
data_for_measurement = malloc(data_size);
- if ( config_noisefilter ) {
+ if (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;
+ 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, 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
@@ -375,40 +390,66 @@ 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 ) {
+ /* 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 ) {
-
+ 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);
-
+ 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;
+ }
+ /* Write Lock */
+ struct flock fl = {F_WRLCK, SEEK_SET, 0, 0, 0};
+ int fd;
+ fl.l_pid = getpid();
+
+ char *outfilename = NULL;
+ chomp(outfile);
+ outfilename = malloc(strlen(outfile) + 1);
+ snprintf(outfilename, LINE_LENGTH - 1, "%s", outfile);
+ 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;
+ 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);
- 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);
+ qargs->n_indexable += pargs->indexable;
+ qargs->n_processed++;
/* Only free cell if found */
cell_free(image.indexed_cell);
@@ -422,54 +463,6 @@ static void process_image(void *pp, int cookie)
}
-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()
@@ -493,34 +486,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)
@@ -554,14 +519,13 @@ 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;
int config_verbose = 0;
int config_satcorr = 1;
int config_checkprefix = 1;
- int config_closer = 1;
+ int config_closer = 0;
int config_insane = 0;
int config_bgsub = 1;
int config_basename = 0;
@@ -585,11 +549,11 @@ 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 nProcesses = 1;
char *prepare_line;
- char prepare_filename[1024];
+ char prepare_filename[LINE_LENGTH];
struct queue_args qargs;
+ struct index_args pargs;
struct beam_params *beam = NULL;
char *element = NULL;
double nominal_photon_energy;
@@ -630,6 +594,7 @@ int main(int argc, char *argv[])
{"threshold", 1, NULL, 't'},
{"no-check-prefix", 0, &config_checkprefix, 0},
{"no-closer-peak", 0, &config_closer, 0},
+ {"closer-peak", 0, &config_closer, 1},
{"insane", 0, &config_insane, 1},
{"image", 1, NULL, 'e'},
{"basename", 0, &config_basename, 1},
@@ -653,7 +618,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 +648,7 @@ int main(int argc, char *argv[])
break;
case 'j' :
- nthreads = atoi(optarg);
+ nProcesses = atoi(optarg);
break;
case 'g' :
@@ -826,7 +791,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 +823,8 @@ int main(int argc, char *argv[])
}
}
- if ( nthreads == 0 ) {
- ERROR("Invalid number of threads.\n");
+ if ( nProcesses == 0 ) {
+ ERROR("Invalid number of processes.\n");
return 1;
}
@@ -976,8 +940,8 @@ int main(int argc, char *argv[])
}
/* Get first filename and use it to set up the indexing */
- prepare_line = malloc(1024*sizeof(char));
- rval = fgets(prepare_line, 1023, fh);
+ prepare_line = malloc(LINE_LENGTH*sizeof(char));
+ rval = fgets(prepare_line, LINE_LENGTH-1, fh);
if ( rval == NULL ) {
ERROR("Failed to get filename to prepare indexing.\n");
return 1;
@@ -989,7 +953,7 @@ int main(int argc, char *argv[])
free(prepare_line);
prepare_line = tmp;
}
- snprintf(prepare_filename, 1023, "%s%s", prefix, prepare_line);
+ snprintf(prepare_filename, LINE_LENGTH-1, "%s%s", prefix, prepare_line);
qargs.use_this_one_instead = prepare_line;
/* Prepare the indexer */
@@ -1027,7 +991,6 @@ int main(int argc, char *argv[])
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;
@@ -1045,11 +1008,185 @@ int main(int argc, char *argv[])
qargs.n_processed = 0;
qargs.n_indexable_last_stats = 0;
qargs.n_processed_last_stats = 0;
+ qargs.updateReader = 0; /* first process updates */
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);
+ /* Read .lst file */
+ register int i;
+ rewind(fh); /* make sure to read from start */
+
+ /* Clear output file content */
+ char *myOutfilename = NULL;
+ chomp(prefix);
+ chomp(outfile);
+ myOutfilename = malloc(strlen(outfile) + 1);
+ snprintf(myOutfilename, LINE_LENGTH - 1, "%s", outfile);
+ FILE *tfh;
+ tfh = fopen(myOutfilename, "a+");
+ if (tfh == NULL) {
+ ERROR("No output filename\n");
+ }
+ fclose(tfh);
+ qargs.static_args.outfile = outfile;
+ int ready_fd;
+ int buff_count;
+ fd_set fdset,tmpset;
+ char buffR[BUFFER], buffW[BUFFER];
+ int fd_pipeIn[nProcesses][2]; /* Process0 In */
+ int fd_pipeOut[nProcesses][2]; /* Process0 Out */
+ unsigned int opts;
+
+ FD_ZERO(&fdset); /* clear the fd_set */
+ /* set pipeIn as non-blocking */
+ for ( i=0; i<nProcesses; i++ ) {
+ opts = fcntl(fd_pipeIn[i][0], F_GETFL);
+ fcntl(fd_pipeIn[i][0], F_SETFL, opts | O_NONBLOCK);
+ }
+
+ /**** PIPING ****/
+ for ( i=0; i<nProcesses; i++ ) {
+ pipe(fd_pipeIn[i]);
+ pipe(fd_pipeOut[i]);
+ }
+
+ int max_fd = 0;
+ for ( i=0; i<nProcesses; i++ ) {
+ FD_SET(fd_pipeIn[i][0], &fdset);
+ if (fd_pipeIn[i][0] > max_fd) { /* find max_fd */
+ max_fd = fd_pipeIn[i][0];
+ }
+ }
+ max_fd = max_fd+1;
+ /* copy file set to tmpset */
+ memcpy((void *) &tmpset,(void *) &fdset, sizeof(fd_set));
+
+ /**** 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 + 1) {
+ exit(0); /* kill */
+ }
+ batchNum = (int) num;
+
+ /**** PLUMBING ****/
+ if (batchNum == qargs.updateReader) {
+ for ( i=0; i<nProcesses; i++ ) {
+ close(fd_pipeIn[i][1]); /* close all write pipes In */
+ close(fd_pipeOut[i][0]); /* close all read pipes Out */
+ }
+ } else {
+ for ( i=0; i<nProcesses; i++ ) {
+ if (i == batchNum - 1) { /* batchNum = 1,2,3 ... */
+ close(fd_pipeIn[i][0]); /* close read pipe In */
+ close(fd_pipeOut[i][1]); /* close write pipe Out */
+ } else {
+ close(fd_pipeIn[i][0]); // close remaining pipes In
+ close(fd_pipeIn[i][1]);
+ close(fd_pipeOut[i][0]); // close remaining pipes Out
+ close(fd_pipeOut[i][1]);
+ }
+ }
+ }
+ /**** INDEXING ****/
+ double tStart, tEnd;
+ tStart = get_monotonic_seconds();
+ int allDone = 0;
+ if (batchNum == qargs.updateReader){
+ char *nextImage = NULL;
+ for ( i=0; i<nProcesses; i++ ) { /* Send out image to all processes*/
+ nextImage = get_pattern(fh);
+ buff_count = sprintf(buffW, "%s",nextImage);
+ write (fd_pipeOut[i][1], buffW, buff_count);
+ }
+ int nFinished = 0;
+ while (!allDone) {
+ /* select from file set for reading */
+ if ((ready_fd = select(max_fd,&fdset,NULL,NULL,NULL)) < 0)
+ perror("select");
+ if (ready_fd > 0) {
+ for ( i=0; i<nProcesses; i++ ) {
+ /* is in file set that raised flag? */
+ if (FD_ISSET(fd_pipeIn[i][0],&fdset)) {
+ /* read from pipe and return number of bytes read */
+ if ((buff_count=read(fd_pipeIn[i][0],&buffR,BUFFER))<0) {
+ perror("read");
+ } else {
+ qargs.n_indexable += atoi(buffR);
+ qargs.n_processed++;
+ /* write to pipe */
+ if ((nextImage = get_pattern(fh)) == NULL){
+ nFinished++; /* no more images */
+ if ( nFinished == nProcesses )
+ allDone = 1; /* EXIT */
+ } else {
+ /* send out image */
+ buff_count = sprintf(buffW, "%s",nextImage);
+ if (write (fd_pipeOut[i][1], buffW, buff_count)<0)
+ perror("write pipe");
+ }
+ }
+ }
+ }
+ }
+ /* file set is modified, so copy original from tmpset */
+ memcpy((void *) &fdset,(void *) &tmpset, sizeof(fd_set));
+
+ /* Update to screen */
+ double tNow = get_monotonic_seconds();
+ if ( tNow >= 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\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_indexable_last_stats = qargs.n_indexable;
+ qargs.n_processed_last_stats = qargs.n_processed;
+ qargs.t_last_stats = tNow;
+ }
+ }
+ /* close my pipes */
+ for ( i=0; i<nProcesses; i++ ) {
+ close(fd_pipeIn[i][0]);
+ close(fd_pipeOut[i][1]);
+ }
+ tEnd = get_monotonic_seconds();
+ printf("Compute Time: %.2fs\n", tEnd - tStart);
+ } else {
+ while(!allDone){
+ /* read from pipe and return number of bytes read */
+ if ((buff_count=read(fd_pipeOut[batchNum-1][0],&buffR,BUFFER))<0) {
+ perror("read1");
+ } else if (buff_count > 0) {
+ /* process image */
+ pargs.filename = buffR;
+ pargs.indexable = 0;
+ process_image(&qargs, &pargs, batchNum);
+ /* request another image */
+ buff_count = sprintf(buffW, "%d\n", pargs.indexable);
+ if(write (fd_pipeIn[batchNum-1][1], buffW, buff_count)<0)
+ perror("write P0");
+ } else if (buff_count == 0) {
+ allDone = 1; /* EXIT */
+ }
+ }
+ /* close my pipes */
+ close(fd_pipeIn[batchNum-1][1]);
+ close(fd_pipeOut[batchNum-1][0]);
+ }
cleanup_indexing(ipriv);
@@ -1065,8 +1202,9 @@ int main(int argc, char *argv[])
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);
-
+ if (batchNum == qargs.updateReader) {
+ STATUS("There were %i images, of which %i could be indexed.\n",
+ qargs.n_processed, qargs.n_indexable);
+ }
return 0;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 03:50:59 +0000