From e132f0a215392b13bf289cad55f2fece6e193625 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 4 Feb 2013 17:42:40 +0100 Subject: Indexing stuff --- src/indexamajig.c | 25 ++++++++++--------------- 1 file changed, 10 insertions(+), 15 deletions(-) (limited to 'src') diff --git a/src/indexamajig.c b/src/indexamajig.c index accfa788..01377de3 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -243,7 +243,6 @@ int main(int argc, char *argv[]) struct index_args iargs; struct beam_params *beam = NULL; char *element = NULL; - double nominal_photon_energy; int stream_flags = STREAM_INTEGRATED; char *hdf5_peak_path = NULL; struct copy_hdf5_field *copyme; @@ -556,7 +555,14 @@ int main(int argc, char *argv[]) } if ( geometry == NULL ) { - ERROR("You need to specify a geometry file with --geometry\n"); + ERROR("You need to provide a geometry file (please read the" + " manual for more details).\n"); + return 1; + } + + if ( beam == NULL ) { + ERROR("You need to provide a beam parameters file (please read" + " the manual for more details).\n"); return 1; } @@ -582,17 +588,6 @@ int main(int argc, char *argv[]) write_stream_header(ofh, argc, argv); - if ( beam != NULL ) { - nominal_photon_energy = beam->photon_energy; - } 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; - } - /* Get first filename and use it to set up the indexing */ prepare_line = malloc(1024); rval = fgets(prepare_line, 1023, fh); @@ -613,8 +608,8 @@ int main(int argc, char *argv[]) /* Prepare the indexer */ if ( indm != NULL ) { - ipriv = prepare_indexing(indm, cell, prepare_filename, det, - nominal_photon_energy); + ipriv = prepare_indexing(indm, cell, prepare_filename, + det, beam); if ( ipriv == NULL ) { ERROR("Failed to prepare indexing.\n"); return 1; -- cgit v1.2.3 From 29cca07716b48f9e433087f5dbb202165b1897e1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 5 Feb 2013 11:36:44 +0100 Subject: WIP on bringing programs up to date --- src/process_hkl.c | 179 +++++++++++++++++++++++++++--------------------------- 1 file changed, 91 insertions(+), 88 deletions(-) (limited to 'src') diff --git a/src/process_hkl.c b/src/process_hkl.c index 8bfbea2b..e9a8daa5 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -3,12 +3,12 @@ * * Assemble and process FEL Bragg intensities * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Lorenzo Galli * * Authors: - * 2009-2012 Thomas White + * 2009-2013 Thomas White * 2011 Andrew Martin * 2012 Lorenzo Galli * @@ -48,6 +48,8 @@ #include "stream.h" #include "reflist.h" #include "image.h" +#include "crystal.h" +#include "thread-pool.h" static void show_help(const char *s) @@ -170,11 +172,11 @@ static double scale_intensities(RefList *reference, RefList *new, } -static int merge_pattern(RefList *model, struct image *new, RefList *reference, - const SymOpList *sym, - double *hist_vals, signed int hist_h, - signed int hist_k, signed int hist_l, int *hist_n, - int config_nopolar) +static int merge_crystal(RefList *model, struct image *image, Crystal *cr, + RefList *reference, const SymOpList *sym, + double *hist_vals, signed int hist_h, + signed int hist_k, signed int hist_l, int *hist_n, + int config_nopolar) { Reflection *refl; RefListIterator *iter; @@ -184,17 +186,18 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, double csx, csy, csz; if ( reference != NULL ) { - scale = scale_intensities(reference, new->reflections, sym); + scale = scale_intensities(reference, + crystal_get_reflections(cr), sym); } else { scale = 1.0; } if ( isnan(scale) ) return 1; - cell_get_reciprocal(new->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); - for ( refl = first_refl(new->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -226,7 +229,7 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, yl = h*asy + k*bsy + l*csy; zl = h*asz + k*bsz + l*csz; - ool = 1.0 / new->lambda; + ool = 1.0 / image->lambda; tt = angle_between(0.0, 0.0, 1.0, xl, yl, zl+ool); phi = atan2(yl, xl); pa = pow(sin(phi)*sin(tt), 2.0); @@ -260,63 +263,72 @@ static int merge_pattern(RefList *model, struct image *new, RefList *reference, } -static void merge_all(FILE *fh, RefList *model, RefList *reference, - int config_startafter, int config_stopafter, - const SymOpList *sym, - int n_total_patterns, - double *hist_vals, signed int hist_h, - signed int hist_k, signed int hist_l, - int *hist_i, int config_nopolar, int min_measurements) +static void display_progress(int n_images, int n_crystals, int n_crystals_used) +{ + if ( !isatty(STDERR_FILENO) ) return; + if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return; + + pthread_mutex_lock(&stderr_lock); + fprintf(stderr, "\r%i images processed, %i crystals, %i crystals used", + n_images, n_crystals, n_crystals_used); + pthread_mutex_unlock(&stderr_lock); + + fflush(stdout); +} + + + +static int merge_all(Stream *st, RefList *model, RefList *reference, + const SymOpList *sym, + double *hist_vals, signed int hist_h, + signed int hist_k, signed int hist_l, + int *hist_i, int config_nopolar, int min_measurements) { int rval; - int n_patterns = 0; - int n_used = 0; + int n_images = 0; + int n_crystals = 0; + int n_crystals_used = 0; Reflection *refl; RefListIterator *iter; - if ( skip_some_files(fh, config_startafter) ) { - ERROR("Failed to skip first %i files.\n", config_startafter); - return; - } - do { struct image image; + int i; image.det = NULL; /* Get data from next chunk */ - rval = read_chunk(fh, &image); + rval = read_chunk(st, &image); if ( rval ) break; - n_patterns++; + n_images++; - if ( (image.reflections != NULL) && (image.indexed_cell) ) { + for ( i=0; i ", hist_h, hist_k, hist_l); @@ -529,40 +526,46 @@ int main(int argc, char *argv[]) } hist_i = 0; - merge_all(fh, model, NULL, config_startafter, config_stopafter, - sym, n_total_patterns, hist_vals, hist_h, hist_k, hist_l, - &hist_i, config_nopolar, min_measurements); - if ( ferror(fh) ) { - ERROR("Stream read error.\n"); + r = merge_all(st, model, NULL, sym, hist_vals, hist_h, hist_k, hist_l, + &hist_i, config_nopolar, min_measurements); + fprintf(stderr, "\n"); + if ( r ) { + ERROR("Error while reading stream.\n"); return 1; } - rewind(fh); if ( config_scale ) { RefList *reference; - STATUS("Extra pass for scaling...\n"); + if ( rewind_stream(st) ) { - reference = copy_reflist(model); + ERROR("Couldn't rewind stream - scaling cannot be " + "performed.\n"); - reflist_free(model); - model = reflist_new(); + } else { - rewind(fh); + int r; - merge_all(fh, model, reference, - config_startafter, config_stopafter, sym, - n_total_patterns, - hist_vals, hist_h, hist_k, hist_l, &hist_i, - config_nopolar, min_measurements); + STATUS("Extra pass for scaling...\n"); - if ( ferror(fh) ) { - ERROR("Stream read error.\n"); - return 1; - } + reference = copy_reflist(model); + + reflist_free(model); + model = reflist_new(); + + r = merge_all(st, model, reference, sym, + hist_vals, hist_h, hist_k, hist_l, &hist_i, + config_nopolar, min_measurements); + fprintf(stderr, "\n"); + if ( r ) { + ERROR("Error while reading stream.\n"); + return 1; + } + + reflist_free(reference); - reflist_free(reference); + } } @@ -579,7 +582,7 @@ int main(int argc, char *argv[]) write_reflist(output, model); - fclose(fh); + close_stream(st); free(sym); reflist_free(model); -- cgit v1.2.3 From 2ce85a95d86e350b785d206405e97d7317672188 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 5 Feb 2013 16:21:40 +0100 Subject: Indexing pipeline - "done"! --- src/im-sandbox.c | 129 ++++++++++++++++++++++++++++++++---------------------- src/im-sandbox.h | 27 ++++-------- src/indexamajig.c | 128 +++++++++++++++-------------------------------------- 3 files changed, 121 insertions(+), 163 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index c8cceb2b..825f2d6a 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -73,14 +73,27 @@ #define STATS_EVERY_N_SECONDS (5) +/* Information about the indexing process for one pattern */ +struct pattern_args +{ + /* "Input" */ + char *filename; + + /* "Output" */ + int n_crystals; +}; + + struct sandbox { pthread_mutex_t lock; - int n_indexable; int n_processed; - int n_indexable_last_stats; + int n_hadcrystals; + int n_crystals; int n_processed_last_stats; + int n_hadcrystals_last_stats; + int n_crystals_last_stats; int t_last_stats; struct index_args *iargs; @@ -164,31 +177,29 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead, static void process_image(const struct index_args *iargs, - struct pattern_args *pargs, FILE *ofh, + struct pattern_args *pargs, Stream *st, int cookie) { float *data_for_measurement; size_t data_size; - 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; int check; struct hdfile *hdfile; struct image image; + int i; image.features = NULL; image.data = NULL; image.flags = NULL; - image.indexed_cell = NULL; image.copyme = iargs->copyme; - image.reflections = NULL; - image.n_saturated = 0; image.id = cookie; image.filename = pargs->filename; image.beam = iargs->beam; image.det = iargs->det; + image.crystals = NULL; + image.n_crystals = 0; hdfile = hdfile_open(image.filename); if ( hdfile == NULL ) return; @@ -291,28 +302,34 @@ static void process_image(const struct index_args *iargs, free(image.data); image.data = data_for_measurement; - /* Calculate orientation matrix (by magic) */ + /* Index the pattern */ + index_pattern(&image, indm, iargs->ipriv); + + /* Default beam parameters */ image.div = image.beam->divergence; image.bw = image.beam->bandwidth; - image.profile_radius = image.beam->profile_radius; - index_pattern(&image, cell, indm, iargs->cellr, - config_verbose, iargs->ipriv, - iargs->config_insane, iargs->tols); + /* Integrate each crystal's diffraction spots */ + for ( i=0; iindexable = 1; + /* Set default crystal parameter(s) */ + crystal_set_profile_radius(image.crystals[i], + image.beam->profile_radius); if ( iargs->integrate_found ) { - image.reflections = select_intersections(&image, - image.indexed_cell); + reflections = select_intersections(&image, + image.crystals[i]); } else { - image.reflections = find_intersections(&image, - image.indexed_cell); + reflections = find_intersections(&image, + image.crystals[i]); } - if (image.reflections != NULL) { + crystal_set_reflections(image.crystals[i], reflections); + + if ( reflections != NULL ) { + integrate_reflections(&image, iargs->config_closer, iargs->config_bgsub, @@ -323,18 +340,15 @@ static void process_image(const struct index_args *iargs, iargs->integrate_saturated); } - } else { - image.reflections = NULL; } - write_chunk(ofh, &image, hdfile, iargs->stream_flags); - fprintf(ofh, "END\n"); - fflush(ofh); + write_chunk(st, &image, hdfile, iargs->include_peaks, + iargs->include_reflections); - /* Only free cell if found */ - cell_free(image.indexed_cell); + for ( i=0; iipriv); + cleanup_indexing(iargs->indm, iargs->ipriv); free(iargs->indm); free(iargs->ipriv); free_detector_geometry(iargs->det); @@ -577,7 +591,7 @@ static void start_worker_process(struct sandbox *sb, int slot) if ( p == 0 ) { - FILE *sfh; + Stream *st; int j; struct sigaction sa; int r; @@ -614,10 +628,10 @@ static void start_worker_process(struct sandbox *sb, int slot) close(filename_pipe[1]); close(result_pipe[0]); - sfh = fdopen(sb->stream_pipe_write[slot], "w"); + st = open_stream_fd_for_write(sb->stream_pipe_write[slot]); run_work(sb->iargs, filename_pipe[0], result_pipe[1], - sfh, slot); - fclose(sfh); + st, slot); + close_stream(st); //close(filename_pipe[0]); close(result_pipe[1]); @@ -697,7 +711,7 @@ static void handle_zombie(struct sandbox *sb) void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, char *use_this_one_instead, - FILE *ofh) + Stream *st) { int i; int allDone; @@ -712,13 +726,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, return; } - sb->n_indexable = 0; sb->n_processed = 0; - sb->n_indexable_last_stats = 0; + sb->n_hadcrystals = 0; + sb->n_crystals = 0; sb->n_processed_last_stats = 0; + sb->n_hadcrystals_last_stats = 0; + sb->n_crystals_last_stats = 0; sb->t_last_stats = get_monotonic_seconds(); sb->n_proc = n_proc; - sb->ofh = ofh; sb->iargs = iargs; pthread_mutex_init(&sb->lock, NULL); @@ -895,8 +910,14 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("Invalid result '%s'\n", results); } } else { - sb->n_indexable += atoi(results); + + int nc = atoi(results); + sb->n_crystals += nc; + if ( nc > 0 ) { + sb->n_hadcrystals++; + } sb->n_processed++; + } /* Send next filename */ @@ -929,14 +950,19 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, tNow = get_monotonic_seconds(); if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) { - STATUS("%i out of %i indexed so far," - " %i out of %i since the last message.\n", - sb->n_indexable, sb->n_processed, - sb->n_indexable - sb->n_indexable_last_stats, - sb->n_processed - sb->n_processed_last_stats); + STATUS("Total so far: %i images processed, " + "%i had crystals, %i crystals overall. " + "Since the last message: %i images processed," + "%i had crystals, %i crystals overall.\n", + sb->n_processed, sb->n_hadcrystals, + sb->n_crystals, + sb->n_processed - sb->n_processed_last_stats, + sb->n_hadcrystals - sb->n_hadcrystals_last_stats, + sb->n_crystals - sb->n_crystals_last_stats); - sb->n_indexable_last_stats = sb->n_indexable; sb->n_processed_last_stats = sb->n_processed; + sb->n_hadcrystals_last_stats = sb->n_hadcrystals; + sb->n_crystals_last_stats = sb->n_crystals; sb->t_last_stats = tNow; } @@ -975,9 +1001,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, free(sb->pids); free(sb->running); - if ( ofh != stdout ) fclose(ofh); - - STATUS("There were %i images, of which %i could be indexed.\n", - sb->n_processed, sb->n_indexable); + STATUS("Final:" + " %i images processed, %i had crystals, %i crystals overall.\n", + sb->n_processed, sb->n_hadcrystals, sb->n_crystals); } diff --git a/src/im-sandbox.h b/src/im-sandbox.h index 872f84ad..176bfde0 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -3,13 +3,13 @@ * * Sandbox for indexing * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon @@ -31,6 +31,8 @@ * */ +#include "stream.h" + enum { PEAK_ZAEF, PEAK_HDF5, @@ -44,10 +46,8 @@ struct index_args int config_cmfilter; int config_noisefilter; int config_verbose; - int stream_flags; /* What goes into the output? */ int config_satcorr; int config_closer; - int config_insane; int config_bgsub; float threshold; float min_gradient; @@ -57,7 +57,6 @@ struct index_args IndexingMethod *indm; IndexingPrivate **ipriv; int peaks; /* Peak detection method */ - int cellr; float tols[4]; struct beam_params *beam; char *element; @@ -70,21 +69,11 @@ struct index_args int use_saturated; int no_revalidate; int integrate_found; + int include_peaks; + int include_reflections; }; - - -/* Information about the indexing process for one pattern */ -struct pattern_args -{ - /* "Input" */ - char *filename; - - /* "Output" */ - int indexable; -}; - extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, - char *use_this_one_instead, FILE *ofh); + char *use_this_one_instead, Stream *st); diff --git a/src/indexamajig.c b/src/indexamajig.c index 01377de3..3076ea86 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -174,37 +174,13 @@ static void show_help(const char *s) } -static int parse_cell_reduction(const char *scellr, int *err, - int *reduction_needs_cell) -{ - *err = 0; - if ( strcmp(scellr, "none") == 0 ) { - *reduction_needs_cell = 0; - return CELLR_NONE; - } else if ( strcmp(scellr, "reduce") == 0) { - *reduction_needs_cell = 1; - return CELLR_REDUCE; - } else if ( strcmp(scellr, "compare") == 0) { - *reduction_needs_cell = 1; - return CELLR_COMPARE; - } else if ( strcmp(scellr, "compare_ab") == 0) { - *reduction_needs_cell = 1; - return CELLR_COMPARE_AB; - } else { - *err = 1; - *reduction_needs_cell = 0; - return CELLR_NONE; - } -} - - int main(int argc, char *argv[]) { int c; char *filename = NULL; char *outfile = NULL; FILE *fh; - FILE *ofh; + Stream *st; char *rval = NULL; int config_noindex = 0; int config_cmfilter = 0; @@ -213,7 +189,6 @@ int main(int argc, char *argv[]) int config_satcorr = 1; int config_checkprefix = 1; int config_closer = 0; - int config_insane = 0; int config_bgsub = 1; int config_basename = 0; float threshold = 800.0; @@ -224,17 +199,13 @@ int main(int argc, char *argv[]) char *geometry = NULL; IndexingMethod *indm; IndexingPrivate **ipriv; - int indexer_needs_cell; - int reduction_needs_cell; char *indm_str = NULL; UnitCell *cell; char *pdb = NULL; char *prefix = NULL; char *speaks = NULL; - char *scellr = NULL; char *toler = NULL; float tols[4] = {5.0, 5.0, 5.0, 1.5}; /* a,b,c,angles (%,%,%,deg) */ - int cellr; int peaks; int n_proc = 1; char *prepare_line; @@ -243,7 +214,6 @@ int main(int argc, char *argv[]) struct index_args iargs; struct beam_params *beam = NULL; char *element = NULL; - int stream_flags = STREAM_INTEGRATED; char *hdf5_peak_path = NULL; struct copy_hdf5_field *copyme; char *intrad = NULL; @@ -255,6 +225,10 @@ int main(int argc, char *argv[]) int no_revalidate = 0; int integrate_found = 0; + /* For ease of upgrading from old method */ + char *scellr = NULL; + int config_insane = 0; + copyme = new_copy_hdf5_field_list(); if ( copyme == NULL ) { ERROR("Couldn't allocate HDF5 field list.\n"); @@ -290,7 +264,6 @@ int main(int argc, char *argv[]) {"peaks", 1, NULL, 2}, {"cell-reduction", 1, NULL, 3}, {"min-gradient", 1, NULL, 4}, - {"record", 1, NULL, 5}, {"cpus", 1, NULL, 6}, {"cpugroup", 1, NULL, 7}, {"cpuoffset", 1, NULL, 8}, @@ -301,6 +274,8 @@ int main(int argc, char *argv[]) {"tolerance", 1, NULL, 13}, {"int-radius", 1, NULL, 14}, + /* FIXME: Add '--no-peaks' and '--no-reflections' */ + {"integrate-saturated",0, &integrate_saturated,1}, {"use-saturated", 0, &use_saturated, 1}, {"no-revalidate", 0, &no_revalidate, 1}, @@ -376,8 +351,7 @@ int main(int argc, char *argv[]) break; case 5 : - stream_flags = parse_stream_flags(optarg); - if ( stream_flags < 0 ) return 1; + ERROR("The option '--record' is no longer used.\n"); break; case 6 : @@ -425,6 +399,12 @@ int main(int argc, char *argv[]) } + if ( scellr != NULL ) { + ERROR("Old-style indexing method specification.\n"); + /* FIXME: Translate to new style */ + return 1; + } + if ( filename == NULL ) { filename = strdup("-"); } @@ -439,17 +419,6 @@ int main(int argc, char *argv[]) } free(filename); - if ( outfile == NULL ) { - ofh = stdout; - } else { - ofh = fopen(outfile, "w"); - 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"); } @@ -482,53 +451,24 @@ int main(int argc, char *argv[]) return 1; } - if ( (indm_str == NULL) || - ((indm_str != NULL) && (strcmp(indm_str, "none") == 0)) ) { - STATUS("Not indexing anything.\n"); - indexer_needs_cell = 0; - reduction_needs_cell = 0; + if ( indm_str == NULL ) { + + STATUS("You didn't specify an indexing method, so I won't try " + " to index anything.\n" + "If that isn't what you wanted, re-run with" + " --indexing=.\n"); indm = NULL; - cellr = CELLR_NONE; + } else { - if ( indm_str == NULL ) { - STATUS("You didn't specify an indexing method, so I " - " won't try to index anything.\n" - "If that isn't what you wanted, re-run with" - " --indexing=.\n"); - indm = NULL; - indexer_needs_cell = 0; - } else { - indm = build_indexer_list(indm_str, - &indexer_needs_cell); - if ( indm == NULL ) { - ERROR("Invalid indexer list '%s'\n", indm_str); - return 1; - } - free(indm_str); - } - reduction_needs_cell = 0; - if ( scellr == NULL ) { - STATUS("You didn't specify a cell reduction method, so" - " I'm going to use 'reduce'.\n"); - cellr = CELLR_REDUCE; - reduction_needs_cell = 1; - } else { - int err; - cellr = parse_cell_reduction(scellr, &err, - &reduction_needs_cell); - if ( err ) { - ERROR("Unrecognised cell reduction '%s'\n", - scellr); - return 1; - } - free(scellr); + indm = build_indexer_list(indm_str); + if ( indm == NULL ) { + ERROR("Invalid indexer list '%s'\n", indm_str); + return 1; } + free(indm_str); } - /* No indexing -> no reduction */ - if ( indm == NULL ) reduction_needs_cell = 0; - if ( toler != NULL ) { int ttt; ttt = sscanf(toler, "%f,%f,%f,%f", @@ -586,7 +526,12 @@ int main(int argc, char *argv[]) cell = NULL; } - write_stream_header(ofh, argc, argv); + st = open_stream_for_write(outfile); + if ( st == NULL ) { + ERROR("Failed to open stream '%s'\n", outfile); + return 1; + } + free(outfile); /* Get first filename and use it to set up the indexing */ prepare_line = malloc(1024); @@ -609,7 +554,7 @@ int main(int argc, char *argv[]) /* Prepare the indexer */ if ( indm != NULL ) { ipriv = prepare_indexing(indm, cell, prepare_filename, - det, beam); + det, beam, tols); if ( ipriv == NULL ) { ERROR("Failed to prepare indexing.\n"); return 1; @@ -627,9 +572,7 @@ int main(int argc, char *argv[]) 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]; @@ -644,7 +587,6 @@ int main(int argc, char *argv[]) iargs.peaks = peaks; iargs.beam = beam; iargs.element = element; - iargs.stream_flags = stream_flags; iargs.hdf5_peak_path = hdf5_peak_path; iargs.copyme = copyme; iargs.ir_inn = ir_inn; @@ -654,9 +596,11 @@ int main(int argc, char *argv[]) iargs.integrate_saturated = integrate_saturated; iargs.no_revalidate = no_revalidate; iargs.integrate_found = integrate_found; + iargs.include_peaks = 1; + iargs.include_reflections = 1; /* FIXME! */ create_sandbox(&iargs, n_proc, prefix, config_basename, fh, - use_this_one_instead, ofh); + use_this_one_instead, st); free(prefix); -- cgit v1.2.3 From 2f06e6a1abf28a22cf78a6c3d89596bec2a6a8c5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 5 Feb 2013 18:05:30 +0100 Subject: Fixes for indexing pipeline --- src/im-sandbox.c | 50 +++++++++++++++++++++++++------------------------- src/im-sandbox.h | 2 +- src/indexamajig.c | 10 ++++++---- 3 files changed, 32 insertions(+), 30 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 825f2d6a..d476c72d 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -305,6 +305,8 @@ static void process_image(const struct index_args *iargs, /* Index the pattern */ index_pattern(&image, indm, iargs->ipriv); + pargs->n_crystals = image.n_crystals; + /* Default beam parameters */ image.div = image.beam->divergence; image.bw = image.beam->bandwidth; @@ -328,22 +330,20 @@ static void process_image(const struct index_args *iargs, crystal_set_reflections(image.crystals[i], reflections); - if ( reflections != NULL ) { - - integrate_reflections(&image, - iargs->config_closer, - iargs->config_bgsub, - iargs->min_int_snr, - iargs->ir_inn, - iargs->ir_mid, - iargs->ir_out, - iargs->integrate_saturated); - } - } - write_chunk(st, &image, hdfile, iargs->include_peaks, - iargs->include_reflections); + /* Integrate all the crystals at once - need all the crystals so that + * overlaps can be detected. */ + integrate_reflections(&image, iargs->config_closer, + iargs->config_bgsub, + iargs->min_int_snr, + iargs->ir_inn, + iargs->ir_mid, + iargs->ir_out, + iargs->integrate_saturated); + + write_chunk(st, &image, hdfile, + iargs->include_peaks, iargs->include_reflections); for ( i=0; ilock, NULL); + sb->ofh = ofh; sb->stream_pipe_read = calloc(n_proc, sizeof(int)); sb->stream_pipe_write = calloc(n_proc, sizeof(int)); if ( sb->stream_pipe_read == NULL ) { @@ -950,15 +953,12 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, tNow = get_monotonic_seconds(); if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) { - STATUS("Total so far: %i images processed, " + STATUS("%i images processed so far, " "%i had crystals, %i crystals overall. " - "Since the last message: %i images processed," - "%i had crystals, %i crystals overall.\n", + "%i images processed since the last message\n", sb->n_processed, sb->n_hadcrystals, sb->n_crystals, - sb->n_processed - sb->n_processed_last_stats, - sb->n_hadcrystals - sb->n_hadcrystals_last_stats, - sb->n_crystals - sb->n_crystals_last_stats); + sb->n_processed - sb->n_processed_last_stats); sb->n_processed_last_stats = sb->n_processed; sb->n_hadcrystals_last_stats = sb->n_hadcrystals; diff --git a/src/im-sandbox.h b/src/im-sandbox.h index 176bfde0..96311056 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -76,4 +76,4 @@ struct index_args extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, - char *use_this_one_instead, Stream *st); + char *use_this_one_instead, FILE *stream); diff --git a/src/indexamajig.c b/src/indexamajig.c index 3076ea86..5d57d336 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -180,7 +180,7 @@ int main(int argc, char *argv[]) char *filename = NULL; char *outfile = NULL; FILE *fh; - Stream *st; + FILE *ofh; char *rval = NULL; int config_noindex = 0; int config_cmfilter = 0; @@ -264,6 +264,7 @@ int main(int argc, char *argv[]) {"peaks", 1, NULL, 2}, {"cell-reduction", 1, NULL, 3}, {"min-gradient", 1, NULL, 4}, + {"record", 1, NULL, 5}, {"cpus", 1, NULL, 6}, {"cpugroup", 1, NULL, 7}, {"cpuoffset", 1, NULL, 8}, @@ -352,6 +353,7 @@ int main(int argc, char *argv[]) case 5 : ERROR("The option '--record' is no longer used.\n"); + /* FIXME: Translate to new style */ break; case 6 : @@ -526,8 +528,8 @@ int main(int argc, char *argv[]) cell = NULL; } - st = open_stream_for_write(outfile); - if ( st == NULL ) { + ofh = fopen(outfile, "w"); + if ( ofh == NULL ) { ERROR("Failed to open stream '%s'\n", outfile); return 1; } @@ -600,7 +602,7 @@ int main(int argc, char *argv[]) iargs.include_reflections = 1; /* FIXME! */ create_sandbox(&iargs, n_proc, prefix, config_basename, fh, - use_this_one_instead, st); + use_this_one_instead, ofh); free(prefix); -- cgit v1.2.3 From 004be7ba8d405c7d04ac6143c783bfec70d296bb Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 5 Feb 2013 23:52:04 +0100 Subject: WIP on updating partialator --- src/hrs-scaling.h | 5 +- src/partialator.c | 203 ++++++++++++++++++++++++++++---------------------- src/post-refinement.h | 3 +- src/scaling-report.h | 10 +-- 4 files changed, 122 insertions(+), 99 deletions(-) (limited to 'src') diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h index 5425b1ff..80940347 100644 --- a/src/hrs-scaling.h +++ b/src/hrs-scaling.h @@ -35,9 +35,10 @@ #endif -#include "image.h" +#include "crystal.h" +#include "reflist.h" -extern RefList *scale_intensities(struct image *images, int n, +extern RefList *scale_intensities(Crystal **crystals, int n, RefList *reference, int n_threads, int noscale); diff --git a/src/partialator.c b/src/partialator.c index c2f6c299..35e25320 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -87,16 +87,16 @@ static void show_help(const char *s) struct refine_args { RefList *full; - struct image *image; + Crystal *crystal; }; struct queue_args { - int n; + int n_started; int n_done; - int n_total_patterns; - struct image *images; + Crystal **crystals; + int n_crystals; struct refine_args task_defaults; }; @@ -104,10 +104,9 @@ struct queue_args static void refine_image(void *task, int id) { struct refine_args *pargs = task; - struct image *image = pargs->image; - image->id = id; + Crystal *cr = pargs->crystal; - pr_refine(image, pargs->full); + pr_refine(cr, pargs->full); } @@ -119,9 +118,9 @@ static void *get_image(void *vqargs) task = malloc(sizeof(struct refine_args)); memcpy(task, &qargs->task_defaults, sizeof(struct refine_args)); - task->image = &qargs->images[qargs->n]; + task->crystal = qargs->crystals[qargs->n_started]; - qargs->n++; + qargs->n_started++; return task; } @@ -133,12 +132,12 @@ static void done_image(void *vqargs, void *task) qargs->n_done++; - progress_bar(qargs->n_done, qargs->n_total_patterns, "Refining"); + progress_bar(qargs->n_done, qargs->n_crystals, "Refining"); free(task); } -static void refine_all(struct image *images, int n_total_patterns, +static void refine_all(Crystal **crystals, int n_crystals, struct detector *det, RefList *full, int nthreads) { @@ -146,19 +145,19 @@ static void refine_all(struct image *images, int n_total_patterns, struct queue_args qargs; task_defaults.full = full; - task_defaults.image = NULL; + task_defaults.crystal = NULL; qargs.task_defaults = task_defaults; - qargs.n = 0; + qargs.n_started = 0; qargs.n_done = 0; - qargs.n_total_patterns = n_total_patterns; - qargs.images = images; + qargs.n_crystals = n_crystals; + qargs.crystals = crystals; /* Don't have threads which are doing nothing */ - if ( n_total_patterns < nthreads ) nthreads = n_total_patterns; + if ( n_crystals < nthreads ) nthreads = n_crystals; run_threads(nthreads, refine_image, get_image, done_image, - &qargs, n_total_patterns, 0, 0, 0); + &qargs, n_crystals, 0, 0, 0); } @@ -201,13 +200,14 @@ static int select_scalable_reflections(RefList *list, RefList *reference) } -static void select_reflections_for_refinement(struct image *images, int n, +static void select_reflections_for_refinement(Crystal **crystals, int n, RefList *full, int have_reference) { int i; for ( i=0; idet = det; - - if ( read_chunk(fh, cur) != 0 ) { - /* Should not happen, because we counted the patterns - * earlier. */ - ERROR("Failed to read chunk from the input stream.\n"); - return 1; + if ( read_chunk(st, cur) != 0 ) { + break; } /* Won't be needing this, if it exists */ image_feature_list_free(cur->features); cur->features = NULL; - - /* "n_usable_patterns" will not be incremented in this case */ - if ( cur->indexed_cell == NULL ) continue; - - /* Fill in initial estimates of stuff */ cur->div = beam->divergence; cur->bw = beam->bandwidth; cur->width = det->max_fs; cur->height = det->max_ss; - cur->osf = 1.0; - cur->profile_radius = beam->profile_radius; - cur->pr_dud = 0; - - /* Muppet proofing */ cur->data = NULL; cur->flags = NULL; cur->beam = NULL; - /* This is the raw list of reflections */ - as = asymmetric_indices(cur->reflections, sym); - reflist_free(cur->reflections); - cur->reflections = as; + n_images++; - update_partialities(cur); + for ( i=0; in_crystals; i++ ) { - nobs += select_scalable_reflections(cur->reflections, - reference); + Crystal *cr; + Crystal **crystals_new; - progress_bar(i, n_total_patterns-1, "Loading pattern data"); - n_usable_patterns++; + crystals_new = realloc(crystals, + (n_crystals+1)*sizeof(Crystal *)); + if ( crystals_new == NULL ) { + ERROR("Failed to allocate memory for crystal " + "list.\n"); + return 1; + } + crystals = crystals_new; + crystals[n_crystals] = crystal_new(); + cr = crystals[n_crystals]; - } - fclose(fh); + /* Fill in initial estimates of stuff */ + crystal_set_osf(cr, 1.0); + crystal_set_profile_radius(cr, beam->profile_radius); + crystal_set_user_flag(cr, 0); + + /* This is the raw list of reflections */ + as = asymmetric_indices(crystal_get_reflections(cr), + sym); + crystal_set_reflections(cr, as); + update_partialities(cur, cr); + + nobs += select_scalable_reflections(as, reference); + + n_crystals++; + + } + + display_progress(n_images, n_crystals); + + } while ( 1 ); + + close_stream(st); /* Make initial estimates */ STATUS("Performing initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); - full = scale_intensities(images, n_usable_patterns, reference, + full = scale_intensities(crystals, n_crystals, reference, nthreads, noscale); - sr = sr_titlepage(images, n_usable_patterns, "scaling-report.pdf", + sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf", infile, cmdline); - sr_iteration(sr, 0, images, n_usable_patterns, full); + sr_iteration(sr, 0, crystals, n_crystals, full); /* Iterate */ for ( i=0; ireflections, - reference); + nobs += select_scalable_reflections(rf, reference); } /* Re-estimate all the full intensities */ reflist_free(full); - full = scale_intensities(images, n_usable_patterns, + full = scale_intensities(crystals, n_crystals, reference, nthreads, noscale); - sr_iteration(sr, i+1, images, n_usable_patterns, full); + sr_iteration(sr, i+1, crystals, n_crystals, full); } sr_finish(sr); n_dud = 0; - for ( i=0; i Date: Wed, 6 Feb 2013 18:59:40 +0100 Subject: Stuff for partialator --- src/hrs-scaling.c | 77 ++++++++++++++++++++++++++------------------------- src/partialator.c | 3 +- src/post-refinement.c | 64 +++++++++++++++++++++++------------------- src/post-refinement.h | 2 +- src/scaling-report.c | 70 +++++++++++++++++++++++++--------------------- 5 files changed, 116 insertions(+), 100 deletions(-) (limited to 'src') diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index f86ae0bb..0498a532 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -57,7 +57,7 @@ struct scale_queue_args { RefList *reference; - struct image *images; + Crystal **crystals; int n_started; double max_shift; }; @@ -65,7 +65,7 @@ struct scale_queue_args struct scale_worker_args { - struct image *image; + Crystal *crystal; double shift; RefList *reference; }; @@ -79,7 +79,7 @@ static void *create_scale_job(void *vqargs) wargs = malloc(sizeof(struct scale_worker_args)); wargs->reference = qargs->reference; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -88,20 +88,21 @@ static void *create_scale_job(void *vqargs) static void run_scale_job(void *vwargs, int cookie) { struct scale_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *reference = wargs->reference; Reflection *refl; RefListIterator *iter; double num = 0.0; double den = 0.0; double corr; + const double osf = crystal_get_osf(cr); - if ( image->pr_dud ) { + if ( crystal_get_user_flag(cr) ) { wargs->shift = 0.0; return; } - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -128,8 +129,8 @@ static void run_scale_job(void *vwargs, int cookie) Ihl = get_intensity(refl) / get_partiality(refl); esd = get_esd_intensity(refl) / get_partiality(refl); - num += Ih * (Ihl/image->osf) / pow(esd/image->osf, 2.0); - den += pow(Ih, 2.0)/pow(esd/image->osf, 2.0); + num += Ih * (Ihl/osf) / pow(esd/osf, 2.0); + den += pow(Ih, 2.0)/pow(esd/osf, 2.0); } @@ -138,7 +139,7 @@ static void run_scale_job(void *vwargs, int cookie) corr = num / den; if ( !isnan(corr) && !isinf(corr) ) { - image->osf *= corr; + crystal_set_osf(cr, osf*corr); } wargs->shift = fabs(corr-1.0); @@ -155,7 +156,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs) } -static double iterate_scale(struct image *images, int n, RefList *reference, +static double iterate_scale(Crystal **crystals, int n, RefList *reference, int n_threads) { struct scale_queue_args qargs; @@ -164,7 +165,7 @@ static double iterate_scale(struct image *images, int n, RefList *reference, qargs.reference = reference; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; qargs.max_shift = 0.0; run_threads(n_threads, run_scale_job, create_scale_job, @@ -178,14 +179,14 @@ struct merge_queue_args { RefList *full; pthread_mutex_t full_lock; - struct image *images; + Crystal **crystals; int n_started; }; struct merge_worker_args { - struct image *image; + Crystal *crystal; RefList *full; pthread_mutex_t *full_lock; }; @@ -200,7 +201,7 @@ static void *create_merge_job(void *vqargs) wargs->full = qargs->full; wargs->full_lock = &qargs->full_lock; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -209,17 +210,17 @@ static void *create_merge_job(void *vqargs) static void run_merge_job(void *vwargs, int cookie) { struct merge_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *full = wargs->full; Reflection *refl; RefListIterator *iter; double G; - if ( image->pr_dud ) return; + if ( crystal_get_user_flag(cr)) return; - G = image->osf; + G = crystal_get_osf(cr); - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -272,7 +273,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs) } -static RefList *lsq_intensities(struct image *images, int n, int n_threads) +static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads) { RefList *full; struct merge_queue_args qargs; @@ -283,7 +284,7 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads) qargs.full = full; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; pthread_mutex_init(&qargs.full_lock, NULL); run_threads(n_threads, run_merge_job, create_merge_job, @@ -308,14 +309,14 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads) struct esd_queue_args { RefList *full; - struct image *images; + Crystal **crystals; int n_started; }; struct esd_worker_args { - struct image *image; + Crystal *crystal; RefList *full; }; @@ -328,7 +329,7 @@ static void *create_esd_job(void *vqargs) wargs = malloc(sizeof(struct esd_worker_args)); wargs->full = qargs->full; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -337,17 +338,17 @@ static void *create_esd_job(void *vqargs) static void run_esd_job(void *vwargs, int cookie) { struct esd_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *full = wargs->full; Reflection *refl; RefListIterator *iter; double G; - if ( image->pr_dud ) return; + if ( crystal_get_user_flag(cr) ) return; - G = image->osf; + G = crystal_get_osf(cr); - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -383,7 +384,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs) } -static void calculate_esds(struct image *images, int n, RefList *full, +static void calculate_esds(Crystal **crystals, int n, RefList *full, int n_threads, int min_red) { struct esd_queue_args qargs; @@ -392,7 +393,7 @@ static void calculate_esds(struct image *images, int n, RefList *full, qargs.full = full; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; for ( refl = first_refl(full, &iter); refl != NULL; @@ -423,7 +424,7 @@ static void calculate_esds(struct image *images, int n, RefList *full, /* Scale the stack of images */ -RefList *scale_intensities(struct image *images, int n, RefList *gref, +RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, int n_threads, int noscale) { int i; @@ -431,17 +432,17 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, RefList *full = NULL; const int min_redundancy = 3; - for ( i=0; i create an initial list to refine against */ if ( gref == NULL ) { - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } /* Iterate */ @@ -457,14 +458,14 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, reference = full; } - max_corr = iterate_scale(images, n, reference, n_threads); + max_corr = iterate_scale(crystals, n, reference, n_threads); //STATUS("Scaling iteration %2i: max correction = %5.2f\n", // i+1, max_corr); /* No reference -> generate list for next iteration */ if ( gref == NULL ) { reflist_free(full); - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } //show_scale_factors(images, n); @@ -474,10 +475,10 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, } while ( (max_corr > 0.01) && (i < MAX_CYCLES) ); if ( gref != NULL ) { - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } /* else we already did it */ - calculate_esds(images, n, full, n_threads, min_redundancy); + calculate_esds(crystals, n, full, n_threads, min_redundancy); return full; } diff --git a/src/partialator.c b/src/partialator.c index 35e25320..30f9bfe4 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -506,6 +506,7 @@ int main(int argc, char *argv[]) crystals = crystals_new; crystals[n_crystals] = crystal_new(); cr = crystals[n_crystals]; + crystal_set_image(cr, cur); /* Fill in initial estimates of stuff */ crystal_set_osf(cr, 1.0); @@ -516,7 +517,7 @@ int main(int argc, char *argv[]) as = asymmetric_indices(crystal_get_reflections(cr), sym); crystal_set_reflections(cr, as); - update_partialities(cur, cr); + update_partialities(cr); nobs += select_scalable_reflections(as, reference); diff --git a/src/post-refinement.c b/src/post-refinement.c index 7410d931..31558ec7 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * * This file is part of CrystFEL. * @@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(struct image *image, int k, Reflection *refl, double r) +double gradient(Crystal *cr, int k, Reflection *refl, double r) { double ds, azix, aziy; double ttlow, tthigh, tt; @@ -103,17 +103,18 @@ double gradient(struct image *image, int k, Reflection *refl, double r) int clamp_low, clamp_high; double klow, khigh; double gr; + struct image *image = crystal_get_image(cr); get_symmetric_indices(refl, &hs, &ks, &ls); - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); xl = hs*asx + ks*bsx + ls*csx; yl = hs*asy + ks*bsy + ls*csy; zl = hs*asz + ks*bsz + ls*csz; - ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls); + ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); @@ -237,8 +238,11 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift) /* Apply the given shift to the 'k'th parameter of 'image'. */ -static void apply_shift(struct image *image, int k, double shift) +static void apply_shift(Crystal *cr, int k, double shift) { + double t; + struct image *image = crystal_get_image(cr); + switch ( k ) { case REF_DIV : @@ -251,7 +255,9 @@ static void apply_shift(struct image *image, int k, double shift) break; case REF_R : - image->profile_radius += shift; + t = crystal_get_profile_radius(cr); + t += shift; + crystal_set_profile_radius(cr, t); break; case REF_ASX : @@ -263,7 +269,7 @@ static void apply_shift(struct image *image, int k, double shift) case REF_CSX : case REF_CSY : case REF_CSZ : - apply_cell_shift(image->indexed_cell, k, shift); + apply_cell_shift(crystal_get_cell(cr), k, shift); break; default : @@ -357,7 +363,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) /* Perform one cycle of post refinement on 'image' against 'full' */ -static double pr_iterate(struct image *image, const RefList *full) +static double pr_iterate(Crystal *cr, const RefList *full) { gsl_matrix *M; gsl_vector *v; @@ -369,7 +375,7 @@ static double pr_iterate(struct image *image, const RefList *full) double max_shift; int nref = 0; - reflections = image->reflections; + reflections = crystal_get_reflections(cr); M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); v = gsl_vector_calloc(NUM_PARAMS); @@ -387,6 +393,7 @@ static double pr_iterate(struct image *image, const RefList *full) double p; Reflection *match; double gradients[NUM_PARAMS]; + const double osf = crystal_get_osf(cr); if ( !get_refinable(refl) ) continue; @@ -398,7 +405,7 @@ static double pr_iterate(struct image *image, const RefList *full) " is it marked as refinable?\n", ha, ka, la); continue; } - I_full = image->osf * get_intensity(match); + I_full = osf * get_intensity(match); /* Actual measurement of this reflection from this pattern? */ I_partial = get_intensity(refl); @@ -412,7 +419,8 @@ static double pr_iterate(struct image *image, const RefList *full) /* Calculate all gradients for this reflection */ for ( k=0; kprofile_radius); + gr = gradient(cr, k, refl, + crystal_get_profile_radius(cr)); gradients[k] = gr; } @@ -429,7 +437,7 @@ static double pr_iterate(struct image *image, const RefList *full) if ( g > k ) continue; M_c = gradients[g] * gradients[k]; - M_c *= w * pow(image->osf * I_full, 2.0); + M_c *= w * pow(osf * I_full, 2.0); M_curr = gsl_matrix_get(M, k, g); gsl_matrix_set(M, k, g, M_curr + M_c); @@ -450,7 +458,7 @@ static double pr_iterate(struct image *image, const RefList *full) //STATUS("%i reflections went into the equations.\n", nref); if ( nref == 0 ) { - image->pr_dud = 1; + crystal_set_user_flag(cr, 1); gsl_matrix_free(M); gsl_vector_free(v); return 0.0; @@ -462,7 +470,7 @@ static double pr_iterate(struct image *image, const RefList *full) for ( param=0; param max_shift ) max_shift = fabs(shift); } @@ -480,7 +488,7 @@ static double pr_iterate(struct image *image, const RefList *full) } -static double guide_dev(struct image *image, const RefList *full) +static double guide_dev(Crystal *cr, const RefList *full) { double dev = 0.0; @@ -488,7 +496,7 @@ static double guide_dev(struct image *image, const RefList *full) Reflection *refl; RefListIterator *iter; - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -508,7 +516,7 @@ static double guide_dev(struct image *image, const RefList *full) * scale_intensities() might not yet have been called, so the * full version may not have been calculated yet. */ - G = image->osf; + G = crystal_get_osf(cr); p = get_partiality(refl); I_partial = get_intensity(refl); I_full = get_intensity(full_version); @@ -524,21 +532,21 @@ static double guide_dev(struct image *image, const RefList *full) } -void pr_refine(struct image *image, const RefList *full) +void pr_refine(Crystal *cr, const RefList *full) { double max_shift, dev; int i; const int verbose = 0; if ( verbose ) { - dev = guide_dev(image, full); + dev = guide_dev(cr, full); STATUS("\n"); /* Deal with progress bar */ STATUS("Before iteration: dev = %10.5e\n", dev); } i = 0; - image->pr_dud = 0; + crystal_set_user_flag(cr, 0); do { double asx, asy, asz; @@ -546,15 +554,15 @@ void pr_refine(struct image *image, const RefList *full) double csx, csy, csz; double dev; - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - max_shift = pr_iterate(image, full); + max_shift = pr_iterate(cr, full); - update_partialities(image); + update_partialities(cr); if ( verbose ) { - dev = guide_dev(image, full); + dev = guide_dev(cr, full); STATUS("PR Iteration %2i: max shift = %10.2f" " dev = %10.5e\n", i+1, max_shift, dev); } diff --git a/src/post-refinement.h b/src/post-refinement.h index 637c0bb0..52e70b2b 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -62,7 +62,7 @@ enum { extern void pr_refine(Crystal *cr, const RefList *full); /* Exported so it can be poked by tests/pr_gradient_check */ -extern double gradient(struct image *image, int k, Reflection *refl, double r); +extern double gradient(Crystal *cr, int k, Reflection *refl, double r); #endif /* POST_REFINEMENT_H */ diff --git a/src/scaling-report.c b/src/scaling-report.c index ba425b96..3de82517 100644 --- a/src/scaling-report.c +++ b/src/scaling-report.c @@ -183,7 +183,7 @@ static void plot_point(cairo_t *cr, double g_width, double g_height, } -static void partiality_graph(cairo_t *cr, const struct image *images, int n, +static void partiality_graph(cairo_t *cr, Crystal **crystals, int n, RefList *full) { const double g_width = 200.0; @@ -219,7 +219,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n, num_nondud = 0; for ( i=0; i osf_max ) osf_max = osf; } osf_max = ceil(osf_max+osf_max/10000.0); @@ -473,9 +477,9 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images, for ( i=0; i= osf_low[b]) @@ -544,8 +548,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images, } -static void intensity_histogram(cairo_t *cr, const struct image *images, - int n, RefList *full, +static void intensity_histogram(cairo_t *cr, Crystal **crystals, int n, + RefList *full, signed int h, signed int k, signed int l) { int f_max; @@ -582,12 +586,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images, for ( i=0; icr); cairo_translate(sr->cr, 480.0, 350.0); - scale_factor_histogram(sr->cr, images, n, + scale_factor_histogram(sr->cr, crystals, n, "Distribution of overall scale factors"); cairo_restore(sr->cr); @@ -830,7 +836,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_save(sr->cr); cairo_translate(sr->cr, 70.0, 330.0); - partiality_graph(sr->cr, images, n, full); + partiality_graph(sr->cr, crystals, n, full); cairo_save(sr->cr); cairo_move_to(sr->cr, 0.0, 0.0); @@ -841,7 +847,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_stroke(sr->cr); cairo_set_dash(sr->cr, NULL, 0, 0.0); cairo_translate(sr->cr, 0.0, -150.0); - partiality_histogram(sr->cr, images, n, full, 1, 0); + partiality_histogram(sr->cr, crystals, n, full, 1, 0); cairo_restore(sr->cr); cairo_save(sr->cr); @@ -854,7 +860,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_set_dash(sr->cr, NULL, 0, 0.0); cairo_translate(sr->cr, 230.0, 200.0); cairo_rotate(sr->cr, -M_PI_2); - partiality_histogram(sr->cr, images, n, full, 0, 1); + partiality_histogram(sr->cr, crystals, n, full, 0, 1); cairo_restore(sr->cr); cairo_restore(sr->cr); @@ -873,7 +879,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_save(sr->cr); cairo_translate(sr->cr, 400.0+140.0*x, 60.0+80.0*y); - intensity_histogram(sr->cr, images, n, full, + intensity_histogram(sr->cr, crystals, n, full, sr->ms_h[i], sr->ms_k[i], sr->ms_l[i]); cairo_restore(sr->cr); -- cgit v1.2.3 From 53382f60c7f8955c016725885c04f510f7fe96ee Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Feb 2013 19:15:06 +0100 Subject: Stuff for partial_sim --- src/partial_sim.c | 66 +++++++++++++++++++++++++++++++++---------------------- 1 file changed, 40 insertions(+), 26 deletions(-) (limited to 'src') diff --git a/src/partial_sim.c b/src/partial_sim.c index facf14ef..4e96fbbc 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -54,11 +54,12 @@ #define NBINS 50 -static void mess_up_cell(UnitCell *cell, double cnoise) +static void mess_up_cell(Crystal *cr, double cnoise) { double ax, ay, az; double bx, by, bz; double cx, cy, cz; + UnitCell *cell = crystal_get_cell(cr); //STATUS("Real:\n"); //cell_print(cell); @@ -82,17 +83,18 @@ static void mess_up_cell(UnitCell *cell, double cnoise) /* For each reflection in "partial", fill in what the intensity would be * according to "full" */ -static void calculate_partials(RefList *partial, double osf, +static void calculate_partials(Crystal *cr, RefList *full, const SymOpList *sym, int random_intensities, pthread_mutex_t *full_lock, unsigned long int *n_ref, double *p_hist, - double *p_max, double max_q, UnitCell *cell) + double *p_max, double max_q) { Reflection *refl; RefListIterator *iter; + double res; - for ( refl = first_refl(partial, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -137,16 +139,17 @@ static void calculate_partials(RefList *partial, double osf, } } - Ip = osf * p * If; + Ip = crystal_get_osf(cr) * p * If; - bin = NBINS*2.0*resolution(cell, h, k, l) / max_q; + res = resolution(crystal_get_cell(cr), h, k, l); + bin = NBINS*2.0*res/max_q; if ( (bin < NBINS) && (bin>=0) ) { p_hist[bin] += p; n_ref[bin]++; if ( p > p_max[bin] ) p_max[bin] = p; } else { STATUS("Reflection out of histogram range: %e %i %f\n", - resolution(cell, h, k, l), bin, p); + res, bin, p); } Ip = gaussian_noise(Ip, 100.0); @@ -206,13 +209,14 @@ struct queue_args unsigned long int n_ref[NBINS]; double p_max[NBINS]; - FILE *stream; + Stream *stream; }; struct worker_args { struct queue_args *qargs; + Crystal *crystal; struct image image; /* Histogram for this image */ @@ -243,21 +247,34 @@ static void *create_job(void *vqargs) static void run_job(void *vwargs, int cookie) { - double osf; struct quaternion orientation; struct worker_args *wargs = vwargs; struct queue_args *qargs = wargs->qargs; int i; + Crystal *cr; + UnitCell *cell; + RefList *reflections; + + cr = crystal_new(); + if ( cr == NULL ) { + ERROR("Failed to create crystal.\n"); + return; + } + wargs->crystal = cr; + crystal_set_image(cr, &wargs->image); - osf = gaussian_noise(1.0, 0.3); + crystal_set_osf(cr, gaussian_noise(1.0, 0.3)); + crystal_set_profile_radius(cr, wargs->image.beam->profile_radius); /* Set up a random orientation */ orientation = random_quaternion(); - wargs->image.indexed_cell = cell_rotate(qargs->cell, orientation); + cell = cell_rotate(qargs->cell, orientation); + crystal_set_cell(cr, cell); + cell_free(cell); snprintf(wargs->image.filename, 255, "dummy.h5"); - wargs->image.reflections = find_intersections(&wargs->image, - wargs->image.indexed_cell); + reflections = find_intersections(&wargs->image, cr); + crystal_set_reflections(cr, reflections); for ( i=0; in_ref[i] = 0; @@ -265,14 +282,14 @@ static void run_job(void *vwargs, int cookie) wargs->p_max[i] = 0.0; } - calculate_partials(wargs->image.reflections, osf, qargs->full, + calculate_partials(cr, qargs->full, qargs->sym, qargs->random_intensities, &qargs->full_lock, wargs->n_ref, wargs->p_hist, wargs->p_max, - qargs->max_q, wargs->image.indexed_cell); + qargs->max_q); /* Give a slightly incorrect cell in the stream */ - mess_up_cell(wargs->image.indexed_cell, qargs->cnoise); + mess_up_cell(cr, qargs->cnoise); } @@ -282,7 +299,7 @@ static void finalise_job(void *vqargs, void *vwargs) struct queue_args *qargs = vqargs; int i; - write_chunk(qargs->stream, &wargs->image, NULL, STREAM_INTEGRATED); + write_chunk(qargs->stream, &wargs->image, NULL, 0, 1); for ( i=0; in_ref[i] += wargs->n_ref[i]; @@ -295,8 +312,7 @@ static void finalise_job(void *vqargs, void *vwargs) qargs->n_done++; progress_bar(qargs->n_done, qargs->n_to_do, "Simulating"); - reflist_free(wargs->image.reflections); - cell_free(wargs->image.indexed_cell); + crystal_free(wargs->crystal); free(wargs); } @@ -315,7 +331,7 @@ int main(int argc, char *argv[]) char *sym_str = NULL; SymOpList *sym; UnitCell *cell = NULL; - FILE *ofh; + Stream *stream; int n = 2; int random_intensities = 0; char *save_file = NULL; @@ -496,13 +512,12 @@ int main(int argc, char *argv[]) ERROR("You must give a filename for the output.\n"); return 1; } - ofh = fopen(output_file, "w"); - if ( ofh == NULL ) { + stream = open_stream_for_write(output_file); + if ( stream == NULL ) { ERROR("Couldn't open output file '%s'\n", output_file); return 1; } free(output_file); - write_stream_header(ofh, argc, argv); image.det = det; image.width = det->max_fs; @@ -511,7 +526,6 @@ int main(int argc, char *argv[]) image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy)); image.div = beam->divergence; image.bw = beam->bandwidth; - image.profile_radius = beam->profile_radius; image.filename = malloc(256); image.copyme = NULL; @@ -528,7 +542,7 @@ int main(int argc, char *argv[]) qargs.random_intensities = random_intensities; qargs.cell = cell; qargs.template_image = ℑ - qargs.stream = ofh; + qargs.stream = stream; qargs.cnoise = cnoise; qargs.max_q = largest_q(&image); @@ -574,7 +588,7 @@ int main(int argc, char *argv[]) } - fclose(ofh); + close_stream(stream); cell_free(cell); free_detector_geometry(det); free(beam); -- cgit v1.2.3 From a1ee07e0887bd23491218301ad440a4ad2efb24f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Feb 2013 20:07:56 +0100 Subject: Stuff --- src/partial_sim.c | 5 +---- src/post-refinement.c | 6 +++--- src/post-refinement.h | 2 +- 3 files changed, 5 insertions(+), 8 deletions(-) (limited to 'src') diff --git a/src/partial_sim.c b/src/partial_sim.c index 4e96fbbc..502e7028 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -252,7 +252,6 @@ static void run_job(void *vwargs, int cookie) struct queue_args *qargs = wargs->qargs; int i; Crystal *cr; - UnitCell *cell; RefList *reflections; cr = crystal_new(); @@ -268,9 +267,7 @@ static void run_job(void *vwargs, int cookie) /* Set up a random orientation */ orientation = random_quaternion(); - cell = cell_rotate(qargs->cell, orientation); - crystal_set_cell(cr, cell); - cell_free(cell); + crystal_set_cell(cr, cell_rotate(qargs->cell, orientation)); snprintf(wargs->image.filename, 255, "dummy.h5"); reflections = find_intersections(&wargs->image, cr); diff --git a/src/post-refinement.c b/src/post-refinement.c index 31558ec7..1439b148 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(Crystal *cr, int k, Reflection *refl, double r) +double gradient(Crystal *cr, int k, Reflection *refl) { double ds, azix, aziy; double ttlow, tthigh, tt; @@ -104,6 +104,7 @@ double gradient(Crystal *cr, int k, Reflection *refl, double r) double klow, khigh; double gr; struct image *image = crystal_get_image(cr); + double r = crystal_get_profile_radius(cr); get_symmetric_indices(refl, &hs, &ks, &ls); @@ -419,8 +420,7 @@ static double pr_iterate(Crystal *cr, const RefList *full) /* Calculate all gradients for this reflection */ for ( k=0; k Date: Wed, 6 Feb 2013 20:39:01 +0100 Subject: Fix process_hkl --- src/partial_sim.c | 5 +++++ src/process_hkl.c | 2 -- 2 files changed, 5 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/partial_sim.c b/src/partial_sim.c index 502e7028..e48c51fa 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -287,6 +287,8 @@ static void run_job(void *vwargs, int cookie) /* Give a slightly incorrect cell in the stream */ mess_up_cell(cr, qargs->cnoise); + + image_add_crystal(&wargs->image, cr); } @@ -523,8 +525,11 @@ int main(int argc, char *argv[]) image.lambda = ph_en_to_lambda(eV_to_J(beam->photon_energy)); image.div = beam->divergence; image.bw = beam->bandwidth; + image.beam = beam; image.filename = malloc(256); image.copyme = NULL; + image.crystals = NULL; + image.n_crystals = 0; if ( random_intensities ) { full = reflist_new(); diff --git a/src/process_hkl.c b/src/process_hkl.c index e9a8daa5..4ba901e6 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -327,8 +327,6 @@ static int merge_all(Stream *st, RefList *model, RefList *reference, } while ( rval == 0 ); - if ( rval ) return 1; - for ( refl = first_refl(model, &iter); refl != NULL; refl = next_refl(refl, iter) ) -- cgit v1.2.3 From df22900f8dba733e8c563704f2686f0e315655f4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Feb 2013 20:56:56 +0100 Subject: Some small(ish) fixes --- src/partialator.c | 4 ++-- src/process_hkl.c | 2 +- src/scaling-report.c | 6 +++--- 3 files changed, 6 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/partialator.c b/src/partialator.c index 30f9bfe4..8b51f826 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -504,7 +504,7 @@ int main(int argc, char *argv[]) return 1; } crystals = crystals_new; - crystals[n_crystals] = crystal_new(); + crystals[n_crystals] = cur->crystals[i]; cr = crystals[n_crystals]; crystal_set_image(cr, cur); @@ -532,7 +532,7 @@ int main(int argc, char *argv[]) close_stream(st); /* Make initial estimates */ - STATUS("Performing initial scaling.\n"); + STATUS("\nPerforming initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); full = scale_intensities(crystals, n_crystals, reference, nthreads, noscale); diff --git a/src/process_hkl.c b/src/process_hkl.c index 4ba901e6..a8abc2e8 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -269,7 +269,7 @@ static void display_progress(int n_images, int n_crystals, int n_crystals_used) if ( tcgetpgrp(STDERR_FILENO) != getpgrp() ) return; pthread_mutex_lock(&stderr_lock); - fprintf(stderr, "\r%i images processed, %i crystals, %i crystals used", + fprintf(stderr, "\r%i images processed, %i crystals, %i crystals used.", n_images, n_crystals, n_crystals_used); pthread_mutex_unlock(&stderr_lock); diff --git a/src/scaling-report.c b/src/scaling-report.c index 3de82517..b161a01e 100644 --- a/src/scaling-report.c +++ b/src/scaling-report.c @@ -3,11 +3,11 @@ * * Write a nice PDF of scaling parameters * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * * This file is part of CrystFEL. * -- cgit v1.2.3 From dac436f21df140b6bf0796f6f9cbb6fcb6c03e2d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Feb 2013 12:35:30 -0800 Subject: Update docs, comments etc --- src/im-sandbox.c | 3 ++- src/indexamajig.c | 53 +++++++++++++++++++++++------------------------------ 2 files changed, 25 insertions(+), 31 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index d476c72d..646b5002 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -954,7 +954,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) { STATUS("%i images processed so far, " - "%i had crystals, %i crystals overall. " + "%i had at least one indexable crystal." + "%i crystals found so far." "%i images processed since the last message\n", sb->n_processed, sb->n_hadcrystals, sb->n_crystals, diff --git a/src/indexamajig.c b/src/indexamajig.c index 5d57d336..8b0f5206 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -3,13 +3,13 @@ * * Index patterns, output hkl+intensity etc. * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon @@ -81,11 +81,7 @@ static void show_help(const char *s) " Default: indexamajig.stream\n" "\n" " --indexing= Use 'methods' for indexing. Provide one or more\n" -" methods separated by commas. Choose from:\n" -" none : no indexing (default)\n" -" dirax : invoke DirAx\n" -" mosflm : invoke MOSFLM (DPS)\n" -" reax : DPS algorithm with known unit cell\n" +" methods separated by commas. See below.\n" " -g. --geometry= Get detector geometry from file.\n" " -b, --beam= Get beam parameters from file (provides nominal\n" " wavelength value if no per-shot value is found in\n" @@ -101,27 +97,7 @@ static void show_help(const char *s) " --hdf5-peaks=

Find peaks table in HDF5 file here.\n" " Default: /processing/hitfinder/peakinfo\n" "\n\n" -"You can control what information is included in the output stream using\n" -"' --record=,,' and so on. Possible flags are:\n\n" -" integrated Include a list of reflection intensities, produced by\n" -" integrating around predicted peak locations.\n" -"\n" -" peaks Include peak locations and intensities from the peak\n" -" search.\n" -"\n" -" peaksifindexed As 'peaks', but only if the pattern could be indexed.\n" -"\n" -" peaksifnotindexed As 'peaks', but only if the pattern could NOT be indexed.\n" -"\n\n" -"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" @@ -157,8 +133,6 @@ static void show_help(const char *s) " --no-check-prefix Don't attempt to correct the --prefix.\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" " integrated intensities.\n" " --use-saturated During the initial peak search, don't reject\n" @@ -170,6 +144,25 @@ static void show_help(const char *s) " --integrate-found Skip the spot prediction step, and just integrate\n" " the intensities of the spots found by the initial\n" " peak search.\n" +"\n" +"Indexing methods:\n\n" +" dirax Invoke DirAx, check linear combinations of the resulting cell\n" +" axes for agreement with your cell, and then check that the cell\n" +" accounts for at least half of the peaks from the peak search.\n" +" mosflm As 'dirax', but invoke MOSFLM instead.\n" +" reax Run the DPS algorithm, looking for the axes of your cell.\n" +"\n" +"You can add the following to the above indexing methods:\n" +" -raw Do not check the resulting unit cell\n" +" (Only for 'dirax' and 'mosflm').\n" +" -axes Check permutations of the axes for correspondence with your cell,\n" +" but do not check linear combinations.\n" +" (Only for 'dirax' and 'mosflm').\n" +" -bad Do not check that the cell accounts for any of the peaks.\n" +"\n" +"The default indexing method is 'none', which means no indexing will be done.\n" +"\n" +"Examples: 'dirax,mosflm,reax', 'dirax-raw,mosflm-raw', 'dirax-raw-bad'\n" ); } -- cgit v1.2.3 From d8ff3bdfcc020ec296d2ad81209fce744343ad57 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 8 Feb 2013 12:36:39 -0800 Subject: Put the command line back in the stream --- src/im-sandbox.c | 10 ++++++---- src/im-sandbox.h | 3 ++- src/indexamajig.c | 2 +- 3 files changed, 9 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 646b5002..f14d5a11 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -569,7 +569,8 @@ static void *run_reader(void *sbv) } -static void start_worker_process(struct sandbox *sb, int slot) +static void start_worker_process(struct sandbox *sb, int slot, + int argc, char *argv[]) { pid_t p; int filename_pipe[2]; @@ -631,6 +632,7 @@ static void start_worker_process(struct sandbox *sb, int slot) close(result_pipe[0]); st = open_stream_fd_for_write(sb->stream_pipe_write[slot]); + write_command(st, argc, argv); run_work(sb->iargs, filename_pipe[0], result_pipe[1], st, slot); close_stream(st); @@ -701,7 +703,7 @@ static void handle_zombie(struct sandbox *sb) STATUS("Last filename was: %s\n", sb->last_filename[i]); sb->n_processed++; - start_worker_process(sb, i); + start_worker_process(sb, i, 0, NULL); } } @@ -713,7 +715,7 @@ static void handle_zombie(struct sandbox *sb) void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, char *use_this_one_instead, - FILE *ofh) + FILE *ofh, int argc, char *argv[]) { int i; int allDone; @@ -823,7 +825,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, /* Fork the right number of times */ lock_sandbox(sb); for ( i=0; i Date: Fri, 8 Feb 2013 15:52:19 -0800 Subject: Small tweaks --- src/im-sandbox.c | 9 ++++----- src/indexamajig.c | 3 ++- 2 files changed, 6 insertions(+), 6 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index f14d5a11..17dc3dc6 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -955,11 +955,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, tNow = get_monotonic_seconds(); if ( tNow >= sb->t_last_stats+STATS_EVERY_N_SECONDS ) { - STATUS("%i images processed so far, " - "%i had at least one indexable crystal." - "%i crystals found so far." - "%i images processed since the last message\n", - sb->n_processed, sb->n_hadcrystals, + STATUS("%4i indexable out of %4i processed, " + "%4i crystals so far. " + "%4i images processed since the last message.\n", + sb->n_hadcrystals, sb->n_processed, sb->n_crystals, sb->n_processed - sb->n_processed_last_stats); diff --git a/src/indexamajig.c b/src/indexamajig.c index d36b2811..97fb1743 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -98,7 +98,7 @@ static void show_help(const char *s) " Default: /processing/hitfinder/peakinfo\n" "\n\n" "For more control over the process, you might need:\n\n" -" --tolerance= Set the tolerances for cell reduction.\n" +" --tolerance= Set the tolerances for cell comparison.\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" @@ -515,6 +515,7 @@ int main(int argc, char *argv[]) return 1; } free(pdb); + STATUS("This is what I understood your unit cell to be:\n"); cell_print(cell); } else { STATUS("No unit cell given.\n"); -- cgit v1.2.3 From afafa4ca344de75f27771d8191e79c6280b68694 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Feb 2013 08:51:43 -0800 Subject: hdfsee: Tidy up confusing error message --- src/dw-hdfsee.c | 1 - src/hdfsee.c | 6 +----- 2 files changed, 1 insertion(+), 6 deletions(-) (limited to 'src') diff --git a/src/dw-hdfsee.c b/src/dw-hdfsee.c index 19aa6f35..ee3b55f0 100644 --- a/src/dw-hdfsee.c +++ b/src/dw-hdfsee.c @@ -1819,7 +1819,6 @@ DisplayWindow *displaywindow_open(const char *filename, const char *peaks, dw->hdfile = hdfile_open(filename); if ( dw->hdfile == NULL ) { - ERROR("Couldn't open file '%s'\n", filename); free(dw); return NULL; } else { diff --git a/src/hdfsee.c b/src/hdfsee.c index be1f4876..894b8af9 100644 --- a/src/hdfsee.c +++ b/src/hdfsee.c @@ -266,11 +266,7 @@ int main(int argc, char *argv[]) ring_radii, n_rings, ring_size); - if ( main_window_list[i] == NULL ) { - ERROR("Couldn't open display window\n"); - } else { - main_n_windows++; - } + if ( main_window_list[i] != NULL ) main_n_windows++; } if ( main_n_windows == 0 ) return 0; -- cgit v1.2.3 From 8907b7cb333a893720cac1def3d86dbe26600fa8 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Feb 2013 10:18:48 -0800 Subject: Clarify lattice type information --- src/indexamajig.c | 22 ++-------------------- 1 file changed, 2 insertions(+), 20 deletions(-) (limited to 'src') diff --git a/src/indexamajig.c b/src/indexamajig.c index 97fb1743..a89ecc5e 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -81,7 +81,8 @@ static void show_help(const char *s) " Default: indexamajig.stream\n" "\n" " --indexing= Use 'methods' for indexing. Provide one or more\n" -" methods separated by commas. See below.\n" +" methods separated by commas.\n" +" See 'man indexamajig' for details.\n" " -g. --geometry= Get detector geometry from file.\n" " -b, --beam= Get beam parameters from file (provides nominal\n" " wavelength value if no per-shot value is found in\n" @@ -144,25 +145,6 @@ static void show_help(const char *s) " --integrate-found Skip the spot prediction step, and just integrate\n" " the intensities of the spots found by the initial\n" " peak search.\n" -"\n" -"Indexing methods:\n\n" -" dirax Invoke DirAx, check linear combinations of the resulting cell\n" -" axes for agreement with your cell, and then check that the cell\n" -" accounts for at least half of the peaks from the peak search.\n" -" mosflm As 'dirax', but invoke MOSFLM instead.\n" -" reax Run the DPS algorithm, looking for the axes of your cell.\n" -"\n" -"You can add the following to the above indexing methods:\n" -" -raw Do not check the resulting unit cell\n" -" (Only for 'dirax' and 'mosflm').\n" -" -axes Check permutations of the axes for correspondence with your cell,\n" -" but do not check linear combinations.\n" -" (Only for 'dirax' and 'mosflm').\n" -" -bad Do not check that the cell accounts for any of the peaks.\n" -"\n" -"The default indexing method is 'none', which means no indexing will be done.\n" -"\n" -"Examples: 'dirax,mosflm,reax', 'dirax-raw,mosflm-raw', 'dirax-raw-bad'\n" ); } -- cgit v1.2.3 From ac00ba6fabcd5fa64b73d79b5a6b783e8763238f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Feb 2013 15:30:10 -0800 Subject: Restore control over what goes in the stream --- src/indexamajig.c | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/indexamajig.c b/src/indexamajig.c index a89ecc5e..57d599b5 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -145,6 +145,8 @@ static void show_help(const char *s) " --integrate-found Skip the spot prediction step, and just integrate\n" " the intensities of the spots found by the initial\n" " peak search.\n" +" --no-peaks-in-stream Do not record peak search results in the stream.\n" +" --no-refls-in-stream Do not record integrated reflections in the stream.\n" ); } @@ -199,6 +201,8 @@ int main(int argc, char *argv[]) int use_saturated = 0; int no_revalidate = 0; int integrate_found = 0; + int config_nopeaks = 0; + int config_norefls = 0; /* For ease of upgrading from old method */ char *scellr = NULL; @@ -249,6 +253,8 @@ int main(int argc, char *argv[]) {"min-integration-snr",1, NULL, 12}, {"tolerance", 1, NULL, 13}, {"int-radius", 1, NULL, 14}, + {"no-peaks-in-stream", 1, &config_nopeaks, 1}, + {"no-refls-in-stream", 1, &config_norefls, 1}, /* FIXME: Add '--no-peaks' and '--no-reflections' */ @@ -574,8 +580,8 @@ int main(int argc, char *argv[]) iargs.integrate_saturated = integrate_saturated; iargs.no_revalidate = no_revalidate; iargs.integrate_found = integrate_found; - iargs.include_peaks = 1; - iargs.include_reflections = 1; /* FIXME! */ + iargs.include_peaks = !config_nopeaks; + iargs.include_reflections = !config_norefls; create_sandbox(&iargs, n_proc, prefix, config_basename, fh, use_this_one_instead, ofh, argc, argv); -- cgit v1.2.3 From 463ab1b598c3ca8a120989f8a26c480a67b66268 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Feb 2013 17:32:27 -0800 Subject: uncenter_cell: Free transformation if it's not needed by the caller --- src/im-sandbox.h | 1 + src/indexamajig.c | 5 +++++ 2 files changed, 6 insertions(+) (limited to 'src') diff --git a/src/im-sandbox.h b/src/im-sandbox.h index 540312ed..80ffc0ad 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -71,6 +71,7 @@ struct index_args int integrate_found; int include_peaks; int include_reflections; + int res_cutoff; }; diff --git a/src/indexamajig.c b/src/indexamajig.c index 57d599b5..908c1b6e 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -121,6 +121,8 @@ static void show_help(const char *s) "-e, --image= Use this image from the HDF5 file.\n" " Example: /data/data0.\n" " Default: The first one found.\n" +" --res-cutoff Estimate the resolution limit of each pattern, and\n" +" don't integrate reflections further out.\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" @@ -203,6 +205,7 @@ int main(int argc, char *argv[]) int integrate_found = 0; int config_nopeaks = 0; int config_norefls = 0; + int config_rescutoff = 0; /* For ease of upgrading from old method */ char *scellr = NULL; @@ -255,6 +258,7 @@ int main(int argc, char *argv[]) {"int-radius", 1, NULL, 14}, {"no-peaks-in-stream", 1, &config_nopeaks, 1}, {"no-refls-in-stream", 1, &config_norefls, 1}, + {"res-cutoff", 1, &config_rescutoff, 1}, /* FIXME: Add '--no-peaks' and '--no-reflections' */ @@ -582,6 +586,7 @@ int main(int argc, char *argv[]) iargs.integrate_found = integrate_found; iargs.include_peaks = !config_nopeaks; iargs.include_reflections = !config_norefls; + iargs.res_cutoff = config_rescutoff; create_sandbox(&iargs, n_proc, prefix, config_basename, fh, use_this_one_instead, ofh, argc, argv); -- cgit v1.2.3 From 5a98496730bbbee2bd5e6da6eff416bd9c8c051f Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Feb 2013 17:35:20 -0800 Subject: Add '--res-cutoff' --- src/im-sandbox.c | 3 ++- 1 file changed, 2 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 17dc3dc6..7494a527 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -340,7 +340,8 @@ static void process_image(const struct index_args *iargs, iargs->ir_inn, iargs->ir_mid, iargs->ir_out, - iargs->integrate_saturated); + iargs->integrate_saturated, + iargs->res_cutoff); write_chunk(st, &image, hdfile, iargs->include_peaks, iargs->include_reflections); -- cgit v1.2.3 From be496cfe14dad210a5b32f387e75b6e2f01ce7a3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Feb 2013 17:36:07 -0800 Subject: Rework indexamajig option processing --- src/im-sandbox.c | 2 +- src/im-sandbox.h | 10 +-- src/indexamajig.c | 226 +++++++++++++++++++++++------------------------------- 3 files changed, 101 insertions(+), 137 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 7494a527..cefdeeb9 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -344,7 +344,7 @@ static void process_image(const struct index_args *iargs, iargs->res_cutoff); write_chunk(st, &image, hdfile, - iargs->include_peaks, iargs->include_reflections); + iargs->stream_peaks, iargs->stream_refls); for ( i=0; i Date: Sat, 9 Feb 2013 17:55:55 -0800 Subject: Fix options --- src/indexamajig.c | 29 ++++++++++++++++------------- 1 file changed, 16 insertions(+), 13 deletions(-) (limited to 'src') diff --git a/src/indexamajig.c b/src/indexamajig.c index 980f23f1..4d0b227a 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -219,28 +219,40 @@ int main(int argc, char *argv[]) /* Long options */ const struct option longopts[] = { + + /* Options with long and short versions */ {"help", 0, NULL, 'h'}, {"input", 1, NULL, 'i'}, {"output", 1, NULL, 'o'}, {"indexing", 1, NULL, 'z'}, {"geometry", 1, NULL, 'g'}, {"beam", 1, NULL, 'b'}, + {"pdb", 1, NULL, 'p'}, + {"prefix", 1, NULL, 'x'}, + {"threshold", 1, NULL, 't'}, + {"image", 1, NULL, 'e'}, + + /* Long-only options with no arguments */ {"filter-cm", 0, &iargs.config_cmfilter, 1}, {"filter-noise", 0, &iargs.config_noisefilter, 1}, {"verbose", 0, &iargs.config_verbose, 1}, - {"pdb", 1, NULL, 'p'}, - {"prefix", 1, NULL, 'x'}, {"no-sat-corr", 0, &iargs.config_satcorr, 0}, {"sat-corr", 0, &iargs.config_satcorr, 1}, - {"threshold", 1, NULL, 't'}, {"no-check-prefix", 0, &config_checkprefix, 0}, {"no-closer-peak", 0, &iargs.config_closer, 0}, {"closer-peak", 0, &iargs.config_closer, 1}, - {"image", 1, NULL, 'e'}, {"basename", 0, &config_basename, 1}, {"bg-sub", 0, &iargs.config_bgsub, 1}, {"no-bg-sub", 0, &iargs.config_bgsub, 0}, + {"no-peaks-in-stream", 0, &iargs.stream_peaks, 0}, + {"no-refls-in-stream", 0, &iargs.stream_refls, 0}, + {"res-cutoff", 0, &iargs.res_cutoff, 1}, + {"integrate-saturated",0, &iargs.integrate_saturated,1}, + {"use-saturated", 0, &iargs.use_saturated, 1}, + {"no-revalidate", 0, &iargs.no_revalidate, 1}, + {"integrate-found", 0, &iargs.integrate_found, 1}, + /* Long-only options with arguments */ {"peaks", 1, NULL, 2}, {"cell-reduction", 1, NULL, 3}, {"min-gradient", 1, NULL, 4}, @@ -254,16 +266,7 @@ int main(int argc, char *argv[]) {"min-integration-snr",1, NULL, 12}, {"tolerance", 1, NULL, 13}, {"int-radius", 1, NULL, 14}, - {"no-peaks-in-stream", 1, &iargs.stream_peaks, 0}, - {"no-refls-in-stream", 1, &iargs.stream_refls, 0}, - {"res-cutoff", 1, &iargs.res_cutoff, 1}, - - /* FIXME: Add '--no-peaks' and '--no-reflections' */ - {"integrate-saturated",0, &iargs.integrate_saturated,1}, - {"use-saturated", 0, &iargs.use_saturated, 1}, - {"no-revalidate", 0, &iargs.no_revalidate, 1}, - {"integrate-found", 0, &iargs.integrate_found, 1}, {0, 0, NULL, 0} }; -- cgit v1.2.3 From 11fe97d53a627fd5caf067bd9fe8f967cf3a72df Mon Sep 17 00:00:00 2001 From: Thomas White Date: Sat, 9 Feb 2013 18:13:15 -0800 Subject: Fussiness, formatting and SPOT --- src/im-sandbox.c | 15 +++++---------- src/im-sandbox.h | 12 ++++++------ src/indexamajig.c | 40 ++++++++++++++++++++-------------------- 3 files changed, 31 insertions(+), 36 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index cefdeeb9..69b22cfc 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -182,9 +182,6 @@ static void process_image(const struct index_args *iargs, { float *data_for_measurement; size_t data_size; - int config_cmfilter = iargs->config_cmfilter; - int config_noisefilter = iargs->config_noisefilter; - IndexingMethod *indm = iargs->indm; int check; struct hdfile *hdfile; struct image image; @@ -258,16 +255,14 @@ static void process_image(const struct index_args *iargs, return; } - if ( config_cmfilter ) { - filter_cm(&image); - } + if ( iargs->cmfilter ) filter_cm(&image); /* Take snapshot of image after CM subtraction but before * the aggressive noise filter. */ data_size = image.width * image.height * sizeof(float); data_for_measurement = malloc(data_size); - if ( config_noisefilter ) { + if ( iargs->noisefilter ) { filter_noise(&image, data_for_measurement); } else { memcpy(data_for_measurement, image.data, data_size); @@ -303,7 +298,7 @@ static void process_image(const struct index_args *iargs, image.data = data_for_measurement; /* Index the pattern */ - index_pattern(&image, indm, iargs->ipriv); + index_pattern(&image, iargs->indm, iargs->ipriv); pargs->n_crystals = image.n_crystals; @@ -334,8 +329,8 @@ static void process_image(const struct index_args *iargs, /* Integrate all the crystals at once - need all the crystals so that * overlaps can be detected. */ - integrate_reflections(&image, iargs->config_closer, - iargs->config_bgsub, + integrate_reflections(&image, iargs->closer, + iargs->bgsub, iargs->min_int_snr, iargs->ir_inn, iargs->ir_mid, diff --git a/src/im-sandbox.h b/src/im-sandbox.h index f7951c04..a454d867 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -43,12 +43,12 @@ enum { struct index_args { UnitCell *cell; - int config_cmfilter; - int config_noisefilter; - int config_verbose; - int config_satcorr; - int config_closer; - int config_bgsub; + int cmfilter; + int noisefilter; + int verbose; + int satcorr; + int closer; + int bgsub; float threshold; float min_gradient; float min_snr; diff --git a/src/indexamajig.c b/src/indexamajig.c index 4d0b227a..ee173c6a 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -179,12 +179,12 @@ int main(int argc, char *argv[]) /* Defaults */ iargs.cell = NULL; - iargs.config_cmfilter = 0; - iargs.config_noisefilter = 0; - iargs.config_verbose = 0; - iargs.config_satcorr = 1; - iargs.config_closer = 0; - iargs.config_bgsub = 1; + iargs.cmfilter = 0; + iargs.noisefilter = 0; + iargs.verbose = 0; + iargs.satcorr = 1; + iargs.closer = 0; + iargs.bgsub = 1; iargs.tols[0] = 5.0; iargs.tols[1] = 5.0; iargs.tols[2] = 5.0; @@ -233,20 +233,20 @@ int main(int argc, char *argv[]) {"image", 1, NULL, 'e'}, /* Long-only options with no arguments */ - {"filter-cm", 0, &iargs.config_cmfilter, 1}, - {"filter-noise", 0, &iargs.config_noisefilter, 1}, - {"verbose", 0, &iargs.config_verbose, 1}, - {"no-sat-corr", 0, &iargs.config_satcorr, 0}, - {"sat-corr", 0, &iargs.config_satcorr, 1}, - {"no-check-prefix", 0, &config_checkprefix, 0}, - {"no-closer-peak", 0, &iargs.config_closer, 0}, - {"closer-peak", 0, &iargs.config_closer, 1}, - {"basename", 0, &config_basename, 1}, - {"bg-sub", 0, &iargs.config_bgsub, 1}, - {"no-bg-sub", 0, &iargs.config_bgsub, 0}, - {"no-peaks-in-stream", 0, &iargs.stream_peaks, 0}, - {"no-refls-in-stream", 0, &iargs.stream_refls, 0}, - {"res-cutoff", 0, &iargs.res_cutoff, 1}, + {"filter-cm", 0, &iargs.cmfilter, 1}, + {"filter-noise", 0, &iargs.noisefilter, 1}, + {"verbose", 0, &iargs.verbose, 1}, + {"no-sat-corr", 0, &iargs.satcorr, 0}, + {"sat-corr", 0, &iargs.satcorr, 1}, + {"no-check-prefix", 0, &config_checkprefix, 0}, + {"no-closer-peak", 0, &iargs.closer, 0}, + {"closer-peak", 0, &iargs.closer, 1}, + {"basename", 0, &config_basename, 1}, + {"bg-sub", 0, &iargs.bgsub, 1}, + {"no-bg-sub", 0, &iargs.bgsub, 0}, + {"no-peaks-in-stream", 0, &iargs.stream_peaks, 0}, + {"no-refls-in-stream", 0, &iargs.stream_refls, 0}, + {"res-cutoff", 0, &iargs.res_cutoff, 1}, {"integrate-saturated",0, &iargs.integrate_saturated,1}, {"use-saturated", 0, &iargs.use_saturated, 1}, {"no-revalidate", 0, &iargs.no_revalidate, 1}, -- cgit v1.2.3 From 9b7274a5b5d8a241277d61e980a8ddf922ebf293 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 11 Feb 2013 20:08:14 -0800 Subject: Improve and fix stream marshalling --- src/im-sandbox.c | 19 ++++++++++++------- 1 file changed, 12 insertions(+), 7 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 69b22cfc..8ca5b122 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -3,13 +3,13 @@ * * Sandbox for indexing * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * Copyright © 2012 Richard Kirian * Copyright © 2012 Lorenzo Galli * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * 2011 Richard Kirian * 2012 Lorenzo Galli * 2012 Chunhong Yoon @@ -481,6 +481,11 @@ static int pump_chunk(FILE *fh, FILE *ofh) } + if ( strcmp(line, "FLUSH\n") == 0 ) { + chunk_finished = 1; + continue; + } + fprintf(ofh, "%s", line); if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) { @@ -629,6 +634,7 @@ static void start_worker_process(struct sandbox *sb, int slot, st = open_stream_fd_for_write(sb->stream_pipe_write[slot]); write_command(st, argc, argv); + write_line(st, "FLUSH"); run_work(sb->iargs, filename_pipe[0], result_pipe[1], st, slot); close_stream(st); @@ -860,10 +866,9 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, r = select(fdmax+1, &fds, NULL, NULL, &tv); if ( r == -1 ) { - if ( errno != EINTR ) { - ERROR("select() failed: %s\n", strerror(errno)); - } - continue; + if ( errno == EINTR ) continue; + ERROR("select() failed: %s\n", strerror(errno)); + break; } if ( r == 0 ) continue; /* No progress this time. Try again */ -- cgit v1.2.3 From 4146f3f9631985b37a8125985dd2c32794d43b59 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 12 Feb 2013 03:02:18 -0800 Subject: Don't lose the last chunks from each worker --- src/im-sandbox.c | 16 ++++++++++++++++ 1 file changed, 16 insertions(+) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 8ca5b122..62f7e667 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -104,6 +104,7 @@ struct sandbox FILE **fhs; int *running; + int *waiting; FILE **result_fhs; int *filename_pipes; int *stream_pipe_read; @@ -416,6 +417,8 @@ static void run_work(const struct index_args *iargs, } + write_line(st, "DONE"); + cleanup_indexing(iargs->indm, iargs->ipriv); free(iargs->indm); free(iargs->ipriv); @@ -486,6 +489,10 @@ static int pump_chunk(FILE *fh, FILE *ofh) continue; } + if ( strcmp(line, "DONE\n") == 0 ) { + return 1; + } + fprintf(ofh, "%s", line); if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) { @@ -554,6 +561,7 @@ static void *run_reader(void *sbv) if ( pump_chunk(sb->fhs[i], sb->ofh) ) { sb->running[i] = 0; + sb->waiting[i] = 0; } } @@ -683,6 +691,7 @@ static void handle_zombie(struct sandbox *sb) int status, p; if ( !sb->running[i] ) continue; + if ( sb->waiting[i] ) continue; p = waitpid(sb->pids[i], &status, WNOHANG); @@ -694,6 +703,7 @@ static void handle_zombie(struct sandbox *sb) if ( p == sb->pids[i] ) { sb->running[i] = 0; + sb->waiting[i] = 1; if ( WIFEXITED(status) ) { continue; @@ -775,6 +785,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->result_fhs = calloc(n_proc, sizeof(FILE *)); sb->pids = calloc(n_proc, sizeof(pid_t)); sb->running = calloc(n_proc, sizeof(int)); + sb->waiting = calloc(n_proc, sizeof(int)); sb->fhs = calloc(sb->n_proc, sizeof(FILE *)); if ( sb->filename_pipes == NULL ) { ERROR("Couldn't allocate memory for pipes.\n"); @@ -792,6 +803,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("Couldn't allocate memory for process flags.\n"); return; } + if ( sb->waiting == NULL ) { + ERROR("Couldn't allocate memory for process flags.\n"); + return; + } sb->last_filename = calloc(n_proc, sizeof(char *)); if ( sb->last_filename == NULL ) { @@ -1003,6 +1018,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, free(sb->result_fhs); free(sb->pids); free(sb->running); + free(sb->waiting); STATUS("Final:" " %i images processed, %i had crystals, %i crystals overall.\n", -- cgit v1.2.3 From 0cda27a3fdfbfb45f3433d9dce541adb79a8762a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 15 Feb 2013 10:47:07 +0100 Subject: Fussiness --- src/im-sandbox.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 62f7e667..8da42d32 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -557,7 +557,9 @@ static void *run_reader(void *sbv) if ( !sb->running[i] ) continue; - if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) continue; + if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) { + continue; + } if ( pump_chunk(sb->fhs[i], sb->ofh) ) { sb->running[i] = 0; -- cgit v1.2.3 From bca26df2b4ad88b56752d62a19b289946f262f94 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 18 Feb 2013 09:25:29 +0100 Subject: Whoops: what was labelled as "gradient" was actually "squared gradient" --- src/indexamajig.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/indexamajig.c b/src/indexamajig.c index ee173c6a..54ae0a5c 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -111,7 +111,7 @@ static void show_help(const char *s) " --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" +" --min-gradient= Minimum squared gradient for Zaefferer peak search.\n" " Default: 100,000.\n" " --min-snr= Minimum signal-to-noise ratio for peaks.\n" " Default: 5.\n" -- cgit v1.2.3 From 672b01c166b8fa015750379a3e86176495cbebf5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 18 Feb 2013 10:54:20 +0100 Subject: indexamajig: Remove --verbose option This hasn't done anything at all since the new indexing subsystem, and hasn't done anything useful for a lot longer. --- src/im-sandbox.h | 1 - src/indexamajig.c | 5 +---- 2 files changed, 1 insertion(+), 5 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.h b/src/im-sandbox.h index a454d867..dfddabe6 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -45,7 +45,6 @@ struct index_args UnitCell *cell; int cmfilter; int noisefilter; - int verbose; int satcorr; int closer; int bgsub; diff --git a/src/indexamajig.c b/src/indexamajig.c index 54ae0a5c..b88469b8 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -128,8 +128,7 @@ static void show_help(const char *s) " --copy-hdf5-field Copy the value of field into the stream. You\n" " can use this option as many times as you need.\n" "\n" -"\nOptions for greater performance or verbosity:\n\n" -" --verbose Be verbose about indexing.\n" +"\nOptions for greater performance:\n\n" " -j Run analyses in parallel. Default 1.\n" "\n" "\nOptions you probably won't need:\n\n" @@ -181,7 +180,6 @@ int main(int argc, char *argv[]) iargs.cell = NULL; iargs.cmfilter = 0; iargs.noisefilter = 0; - iargs.verbose = 0; iargs.satcorr = 1; iargs.closer = 0; iargs.bgsub = 1; @@ -235,7 +233,6 @@ int main(int argc, char *argv[]) /* Long-only options with no arguments */ {"filter-cm", 0, &iargs.cmfilter, 1}, {"filter-noise", 0, &iargs.noisefilter, 1}, - {"verbose", 0, &iargs.verbose, 1}, {"no-sat-corr", 0, &iargs.satcorr, 0}, {"sat-corr", 0, &iargs.satcorr, 1}, {"no-check-prefix", 0, &config_checkprefix, 0}, -- cgit v1.2.3 From e7866b44c0704bc76a2bed2c0f8950797a6055f5 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 18 Feb 2013 11:35:36 +0100 Subject: Use a different folder for each worker process --- src/im-sandbox.c | 40 +++++++++++++++++++++++++++++++++++++++- 1 file changed, 39 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 8da42d32..aeeb023c 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -47,6 +47,7 @@ #include #include #include +#include #ifdef HAVE_CLOCK_GETTIME #include @@ -187,13 +188,17 @@ static void process_image(const struct index_args *iargs, struct hdfile *hdfile; struct image image; int i; + char filename[1024]; + + /* Prefix to jump out of temporary folder */ + snprintf(filename, 1023, "../../%s", pargs->filename); image.features = NULL; image.data = NULL; image.flags = NULL; image.copyme = iargs->copyme; image.id = cookie; - image.filename = pargs->filename; + image.filename = filename; image.beam = iargs->beam; image.det = iargs->det; image.crystals = NULL; @@ -580,6 +585,35 @@ static void *run_reader(void *sbv) } +static int create_temporary_folder(signed int n) +{ + int r; + char tmp[64]; + + if ( n < 0 ) { + snprintf(tmp, 63, "indexamajig.%i", getpid()); + } else { + snprintf(tmp, 63, "worker.%i", n); + } + + r = mkdir(tmp, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); + if ( r ) { + ERROR("Failed to create temporary folder: %s\n", + strerror(errno)); + return 1; + } + + r = chdir(tmp); + if ( r ) { + ERROR("Failed to chdir to temporary folder: %s\n", + strerror(errno)); + return 1; + } + + return 0; +} + + static void start_worker_process(struct sandbox *sb, int slot, int argc, char *argv[]) { @@ -620,6 +654,8 @@ static void start_worker_process(struct sandbox *sb, int slot, return; } + create_temporary_folder(slot); + /* Free resources which will not be needed by worker */ for ( j=0; jn_proc; j++ ) { if ( (j != slot) && (sb->running[j]) ) { @@ -841,6 +877,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, return; } + if ( create_temporary_folder(-1) ) return; + /* Fork the right number of times */ lock_sandbox(sb); for ( i=0; i Date: Tue, 19 Feb 2013 12:21:17 +0100 Subject: Don't add ../../ to the filename if it's already absolute --- src/im-sandbox.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index aeeb023c..5d42d17d 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -191,7 +191,9 @@ static void process_image(const struct index_args *iargs, char filename[1024]; /* Prefix to jump out of temporary folder */ - snprintf(filename, 1023, "../../%s", pargs->filename); + if ( pargs->filename[0] != '/' ) { + snprintf(filename, 1023, "../../%s", pargs->filename); + } image.features = NULL; image.data = NULL; -- cgit v1.2.3 From a5439f150d97b30576ea9f337c004ff0a9e404a1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 19 Feb 2013 13:51:47 +0100 Subject: D'oh --- src/im-sandbox.c | 2 ++ 1 file changed, 2 insertions(+) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 5d42d17d..d393c356 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -193,6 +193,8 @@ static void process_image(const struct index_args *iargs, /* Prefix to jump out of temporary folder */ if ( pargs->filename[0] != '/' ) { snprintf(filename, 1023, "../../%s", pargs->filename); + } else { + snprintf(filename, 1023, "%s", pargs->filename); } image.features = NULL; -- cgit v1.2.3 From 8d28f13e79dbec60a6f95a85ecfcf16ba47a73c7 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 19 Feb 2013 18:13:25 +0100 Subject: Fix incorrect path in stream (which broke check-near-bragg et al.) --- src/im-sandbox.c | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index d393c356..0afe9ff3 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -202,13 +202,13 @@ static void process_image(const struct index_args *iargs, image.flags = NULL; image.copyme = iargs->copyme; image.id = cookie; - image.filename = filename; + image.filename = pargs->filename; /* Relative to top level */ image.beam = iargs->beam; image.det = iargs->det; image.crystals = NULL; image.n_crystals = 0; - hdfile = hdfile_open(image.filename); + hdfile = hdfile_open(filename); /* Relative to temporary folder */ if ( hdfile == NULL ) return; if ( iargs->element != NULL ) { -- cgit v1.2.3 From b821aae2dd29cbdec5ccd33a3929a61ee047f0e6 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 10:15:23 +0100 Subject: More robust stream marshalling --- src/im-sandbox.c | 98 ++++++++++++++++++++++---------------------------------- 1 file changed, 39 insertions(+), 59 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 0afe9ff3..37ccc35c 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -105,7 +105,6 @@ struct sandbox FILE **fhs; int *running; - int *waiting; FILE **result_fhs; int *filename_pipes; int *stream_pipe_read; @@ -426,8 +425,6 @@ static void run_work(const struct index_args *iargs, } - write_line(st, "DONE"); - cleanup_indexing(iargs->indm, iargs->ipriv); free(iargs->indm); free(iargs->ipriv); @@ -468,7 +465,6 @@ static time_t get_monotonic_seconds() static int pump_chunk(FILE *fh, FILE *ofh) { int chunk_started = 0; - int chunk_finished = 0; do { @@ -479,40 +475,29 @@ static int pump_chunk(FILE *fh, FILE *ofh) if ( rval == NULL ) { if ( feof(fh) ) { - /* Process died */ + /* Whoops, connection lost */ if ( chunk_started ) { ERROR("EOF during chunk!\n"); fprintf(ofh, "Chunk is unfinished!\n"); - } + fprintf(ofh, CHUNK_END_MARKER"\n"); + } /* else normal end of output */ return 1; } else { ERROR("fgets() failed: %s\n", strerror(errno)); } - chunk_finished = 1; - continue; + break; } - if ( strcmp(line, "FLUSH\n") == 0 ) { - chunk_finished = 1; - continue; - } - - if ( strcmp(line, "DONE\n") == 0 ) { - return 1; - } + if ( strcmp(line, "FLUSH\n") == 0 ) break; fprintf(ofh, "%s", line); - if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) { - chunk_finished = 1; - } - if ( strcmp(line, CHUNK_START_MARKER"\n") == 0 ) { - chunk_started = 1; - } + if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) break; + if ( strcmp(line, CHUNK_START_MARKER"\n") == 0 ) break; - } while ( !chunk_finished ); + } while ( 1 ); return 0; } @@ -529,7 +514,7 @@ static void *run_reader(void *sbv) fd_set fds; int fdmax; - tv.tv_sec = 5; + tv.tv_sec = 1; tv.tv_usec = 0; FD_ZERO(&fds); @@ -539,7 +524,9 @@ static void *run_reader(void *sbv) int fd; - if ( !sb->running[i] ) continue; + /* Listen for output from all processes which have a + * connection, even if they're not "running". */ + if ( sb->fhs[i] == NULL ) continue; fd = sb->stream_pipe_read[i]; @@ -559,28 +546,24 @@ static void *run_reader(void *sbv) continue; } - if ( r == 0 ) continue; /* Nothing this time. Try again */ - lock_sandbox(sb); for ( i=0; in_proc; i++ ) { - if ( !sb->running[i] ) continue; - if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) { continue; } + /* If the chunk cannot be read, assume the connection + * is broken and that the process will die soon. */ if ( pump_chunk(sb->fhs[i], sb->ofh) ) { - sb->running[i] = 0; - sb->waiting[i] = 0; + sb->fhs[i] = NULL; } } done = 1; - for ( i=0; in_proc; i++ ) { - if ( sb->running[i] ) done = 0; - } + if ( sb->running != NULL ) done = 0; + unlock_sandbox(sb); } @@ -733,7 +716,6 @@ static void handle_zombie(struct sandbox *sb) int status, p; if ( !sb->running[i] ) continue; - if ( sb->waiting[i] ) continue; p = waitpid(sb->pids[i], &status, WNOHANG); @@ -745,7 +727,6 @@ static void handle_zombie(struct sandbox *sb) if ( p == sb->pids[i] ) { sb->running[i] = 0; - sb->waiting[i] = 1; if ( WIFEXITED(status) ) { continue; @@ -827,7 +808,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->result_fhs = calloc(n_proc, sizeof(FILE *)); sb->pids = calloc(n_proc, sizeof(pid_t)); sb->running = calloc(n_proc, sizeof(int)); - sb->waiting = calloc(n_proc, sizeof(int)); sb->fhs = calloc(sb->n_proc, sizeof(FILE *)); if ( sb->filename_pipes == NULL ) { ERROR("Couldn't allocate memory for pipes.\n"); @@ -845,10 +825,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("Couldn't allocate memory for process flags.\n"); return; } - if ( sb->waiting == NULL ) { - ERROR("Couldn't allocate memory for process flags.\n"); - return; - } sb->last_filename = calloc(n_proc, sizeof(char *)); if ( sb->last_filename == NULL ) { @@ -861,11 +837,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, } unlock_sandbox(sb); - if ( pthread_create(&reader_thread, NULL, run_reader, (void *)sb) ) { - ERROR("Failed to create reader thread.\n"); - return; - } - if ( pipe(signal_pipe) == -1 ) { ERROR("Failed to create signal pipe.\n"); return; @@ -890,6 +861,13 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, } unlock_sandbox(sb); + /* Start reader thread after forking, so that things are definitely + * "running" */ + if ( pthread_create(&reader_thread, NULL, run_reader, (void *)sb) ) { + ERROR("Failed to create reader thread.\n"); + return; + } + allDone = 0; while ( !allDone ) { @@ -909,9 +887,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int fd; - if ( !sb->running[i] ) { - continue; - } + if ( sb->result_fhs[i] == NULL) continue; fd = fileno(sb->result_fhs[i]); FD_SET(fd, &fds); @@ -949,14 +925,10 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int fd; char *eptr; - if ( !sb->running[i] ) { - continue; - } + if ( sb->result_fhs[i] == NULL ) continue; fd = fileno(sb->result_fhs[i]); - if ( !FD_ISSET(fd, &fds) ) { - continue; - } + if ( !FD_ISSET(fd, &fds) ) continue; rval = fgets(results, 1024, sb->result_fhs[i]); if ( rval == NULL ) { @@ -964,6 +936,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("fgets() failed: %s\n", strerror(errno)); } + sb->result_fhs[i] = NULL; continue; } @@ -1042,7 +1015,13 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, fclose(fh); - pthread_mutex_destroy(&sb->lock); + /* Indicate to the reader thread that we are done */ + lock_sandbox(sb); + free(sb->running); + sb->running = NULL; + unlock_sandbox(sb); + + pthread_join(reader_thread, NULL); for ( i=0; ifilename_pipes[i]); - fclose(sb->result_fhs[i]); + if ( sb->result_fhs[i] != NULL ) fclose(sb->result_fhs[i]); } for ( i=0; in_proc; i++ ) { @@ -1061,11 +1040,12 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, free(sb->filename_pipes); free(sb->result_fhs); free(sb->pids); - free(sb->running); - free(sb->waiting); + + pthread_mutex_destroy(&sb->lock); STATUS("Final:" " %i images processed, %i had crystals, %i crystals overall.\n", sb->n_processed, sb->n_hadcrystals, sb->n_crystals); + free(sb); } -- cgit v1.2.3 From e995d87d18a52ff2b18209459b33d2988dd7aab3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 11:40:46 +0100 Subject: Even more robust stream handling --- src/im-sandbox.c | 192 ++++++++++++++++++++++++++++++++++++------------------- 1 file changed, 125 insertions(+), 67 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 37ccc35c..b8f45247 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -85,6 +85,24 @@ struct pattern_args }; +struct sb_reader +{ + pthread_mutex_t lock; + int done; + + /* If a worker process dies unexpectedly (e.g. if it segfaults), then + * the pipe for its output can still stay open for a little while while + * its buffer empties. The number of pipes being read from is therefore + * not necessarily the same as the number of worker processes. */ + int n_read; + FILE **fhs; + int *fds; + + /* Final output file handle */ + FILE *ofh; +}; + + struct sandbox { pthread_mutex_t lock; @@ -101,15 +119,14 @@ struct sandbox int n_proc; pid_t *pids; - FILE *ofh; - FILE **fhs; int *running; FILE **result_fhs; int *filename_pipes; - int *stream_pipe_read; int *stream_pipe_write; char **last_filename; + + struct sb_reader *reader; }; @@ -502,40 +519,96 @@ static int pump_chunk(FILE *fh, FILE *ofh) } -static void *run_reader(void *sbv) +/* Add an fd to the list of pipes to be read from */ +static void add_pipe(struct sb_reader *rd, int fd) { - struct sandbox *sb = sbv; - int done = 0; + int *fds_new; + FILE **fhs_new; + int slot; + + pthread_mutex_lock(&rd->lock); + + fds_new = realloc(rd->fds, (rd->n_read+1)*sizeof(int)); + if ( fds_new == NULL ) { + ERROR("Failed to allocate memory for new pipe.\n"); + return; + } - while ( !done ) { + fhs_new = realloc(rd->fhs, (rd->n_read+1)*sizeof(FILE *)); + if ( fhs_new == NULL ) { + ERROR("Failed to allocate memory for new FH.\n"); + return; + } + + rd->fds = fds_new; + rd->fhs = fhs_new; + slot = rd->n_read; + + rd->fds[slot] = fd; + + rd->fhs[slot] = fdopen(fd, "r"); + if ( rd->fhs[slot] == NULL ) { + ERROR("Couldn't fdopen() stream!\n"); + return; + } + + rd->n_read++; + + pthread_mutex_unlock(&rd->lock); +} + + +/* Assumes that the caller is already holding rd->lock! */ +static void remove_pipe(struct sb_reader *rd, int d) +{ + int i; + + for ( i=d; in_read; i++ ) { + if ( i < rd->n_read-1 ) { + rd->fds[i] = rd->fds[i+1]; + rd->fhs[i] = rd->fhs[i+1]; + } /* else don't bother */ + } + + rd->n_read--; + + /* We don't bother shrinking the arrays */ +} + + +static void *run_reader(void *rdv) +{ + struct sb_reader *rd = rdv; + + while ( 1 ) { int r, i; struct timeval tv; fd_set fds; int fdmax; + /* Exit when: + * - No fhs left open to read from + * AND - Main thread says "done" */ + if ( (rd->n_read == 0) && rd->done ) break; + tv.tv_sec = 1; tv.tv_usec = 0; FD_ZERO(&fds); fdmax = 0; - lock_sandbox(sb); - for ( i=0; in_proc; i++ ) { + pthread_mutex_lock(&rd->lock); + for ( i=0; in_read; i++ ) { int fd; - /* Listen for output from all processes which have a - * connection, even if they're not "running". */ - if ( sb->fhs[i] == NULL ) continue; - - fd = sb->stream_pipe_read[i]; + fd = rd->fds[i]; FD_SET(fd, &fds); if ( fd > fdmax ) fdmax = fd; } - - unlock_sandbox(sb); + pthread_mutex_unlock(&rd->lock); r = select(fdmax+1, &fds, NULL, NULL, &tv); @@ -546,25 +619,23 @@ static void *run_reader(void *sbv) continue; } - lock_sandbox(sb); - for ( i=0; in_proc; i++ ) { + pthread_mutex_lock(&rd->lock); + for ( i=0; in_read; i++ ) { - if ( !FD_ISSET(sb->stream_pipe_read[i], &fds) ) { + if ( !FD_ISSET(rd->fds[i], &fds) ) { continue; } /* If the chunk cannot be read, assume the connection * is broken and that the process will die soon. */ - if ( pump_chunk(sb->fhs[i], sb->ofh) ) { - sb->fhs[i] = NULL; + if ( pump_chunk(rd->fhs[i], rd->ofh) ) { + /* remove_pipe() assumes that the caller is + * holding rd->lock ! */ + remove_pipe(rd, i); } } - - done = 1; - if ( sb->running != NULL ) done = 0; - - unlock_sandbox(sb); + pthread_mutex_unlock(&rd->lock); } @@ -607,6 +678,7 @@ static void start_worker_process(struct sandbox *sb, int slot, pid_t p; int filename_pipe[2]; int result_pipe[2]; + int stream_pipe[2]; if ( pipe(filename_pipe) == - 1 ) { ERROR("pipe() failed!\n"); @@ -618,6 +690,11 @@ static void start_worker_process(struct sandbox *sb, int slot, return; } + if ( pipe(stream_pipe) == - 1 ) { + ERROR("pipe() failed!\n"); + return; + } + p = fork(); if ( p == -1 ) { ERROR("fork() failed!\n"); @@ -665,7 +742,7 @@ static void start_worker_process(struct sandbox *sb, int slot, close(filename_pipe[1]); close(result_pipe[0]); - st = open_stream_fd_for_write(sb->stream_pipe_write[slot]); + st = open_stream_fd_for_write(stream_pipe[1]); write_command(st, argc, argv); write_line(st, "FLUSH"); run_work(sb->iargs, filename_pipe[0], result_pipe[1], @@ -683,14 +760,11 @@ static void start_worker_process(struct sandbox *sb, int slot, * and the 'read' end of the result pipe. */ sb->pids[slot] = p; sb->running[slot] = 1; + add_pipe(sb->reader, stream_pipe[0]); close(filename_pipe[0]); close(result_pipe[1]); + close(stream_pipe[1]); sb->filename_pipes[slot] = filename_pipe[1]; - sb->fhs[slot] = fdopen(sb->stream_pipe_read[slot], "r"); - if ( sb->fhs[slot] == NULL ) { - ERROR("Couldn't fdopen() stream!\n"); - return; - } sb->result_fhs[slot] = fdopen(result_pipe[0], "r"); if ( sb->result_fhs[slot] == NULL ) { @@ -765,6 +839,16 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, return; } + sb->reader = calloc(1, sizeof(struct sb_reader)); + if ( sb->reader == NULL ) { + ERROR("Couldn't allocate memory for SB reader.\n"); + free(sb); + return; + } + + pthread_mutex_init(&sb->lock, NULL); + pthread_mutex_init(&sb->reader->lock, NULL); + sb->n_processed = 0; sb->n_hadcrystals = 0; sb->n_crystals = 0; @@ -775,40 +859,21 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->n_proc = n_proc; sb->iargs = iargs; - pthread_mutex_init(&sb->lock, NULL); + sb->reader->fds = NULL; + sb->reader->fhs = NULL; + sb->reader->ofh = ofh; - sb->ofh = ofh; - sb->stream_pipe_read = calloc(n_proc, sizeof(int)); sb->stream_pipe_write = calloc(n_proc, sizeof(int)); - if ( sb->stream_pipe_read == NULL ) { - ERROR("Couldn't allocate memory for pipes.\n"); - return; - } if ( sb->stream_pipe_write == NULL ) { ERROR("Couldn't allocate memory for pipes.\n"); return; } - for ( i=0; istream_pipe_read[i] = stream_pipe[0]; - sb->stream_pipe_write[i] = stream_pipe[1]; - - } - lock_sandbox(sb); sb->filename_pipes = calloc(n_proc, sizeof(int)); sb->result_fhs = calloc(n_proc, sizeof(FILE *)); sb->pids = calloc(n_proc, sizeof(pid_t)); sb->running = calloc(n_proc, sizeof(int)); - sb->fhs = calloc(sb->n_proc, sizeof(FILE *)); if ( sb->filename_pipes == NULL ) { ERROR("Couldn't allocate memory for pipes.\n"); return; @@ -831,10 +896,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, ERROR("Couldn't allocate memory for last filename list.\n"); return; } - if ( sb->fhs == NULL ) { - ERROR("Couldn't allocate memory for file handles!\n"); - return; - } unlock_sandbox(sb); if ( pipe(signal_pipe) == -1 ) { @@ -863,7 +924,8 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, /* Start reader thread after forking, so that things are definitely * "running" */ - if ( pthread_create(&reader_thread, NULL, run_reader, (void *)sb) ) { + if ( pthread_create(&reader_thread, NULL, run_reader, + (void *)sb->reader) ) { ERROR("Failed to create reader thread.\n"); return; } @@ -1016,10 +1078,9 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, fclose(fh); /* Indicate to the reader thread that we are done */ - lock_sandbox(sb); - free(sb->running); - sb->running = NULL; - unlock_sandbox(sb); + pthread_mutex_lock(&sb->reader->lock); + sb->reader->done = 1; + pthread_mutex_unlock(&sb->reader->lock); pthread_join(reader_thread, NULL); @@ -1033,10 +1094,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, if ( sb->result_fhs[i] != NULL ) fclose(sb->result_fhs[i]); } - for ( i=0; in_proc; i++ ) { - fclose(sb->fhs[i]); - } - free(sb->fhs); + free(sb->running); free(sb->filename_pipes); free(sb->result_fhs); free(sb->pids); -- cgit v1.2.3 From 9dd569008c6190b697d64c77a96e50a59e65fbf0 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 11:40:56 +0100 Subject: Create temporary folders only when necessary They might already exist if a previous worker process died and was restarted. --- src/im-sandbox.c | 19 ++++++++++++++----- 1 file changed, 14 insertions(+), 5 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index b8f45247..504c52c4 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -647,6 +647,7 @@ static int create_temporary_folder(signed int n) { int r; char tmp[64]; + struct stat s; if ( n < 0 ) { snprintf(tmp, 63, "indexamajig.%i", getpid()); @@ -654,11 +655,19 @@ static int create_temporary_folder(signed int n) snprintf(tmp, 63, "worker.%i", n); } - r = mkdir(tmp, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); - if ( r ) { - ERROR("Failed to create temporary folder: %s\n", - strerror(errno)); - return 1; + if ( stat(tmp, &s) == -1 ) { + if ( errno != ENOENT ) { + ERROR("Failed to stat temporary folder.\n"); + return 1; + } + + r = mkdir(tmp, S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); + if ( r ) { + ERROR("Failed to create temporary folder: %s\n", + strerror(errno)); + return 1; + } + } r = chdir(tmp); -- cgit v1.2.3 From 8a1898f66575495b9dbb75a7ce1bc6ea26c05cb3 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 11:57:26 +0100 Subject: Move processing pipeline to separate file --- src/im-sandbox.c | 210 +---------------------------------------------- src/im-sandbox.h | 44 +--------- src/process_image.c | 231 ++++++++++++++++++++++++++++++++++++++++++++++++++++ src/process_image.h | 94 +++++++++++++++++++++ 4 files changed, 329 insertions(+), 250 deletions(-) create mode 100644 src/process_image.c create mode 100644 src/process_image.h (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 504c52c4..290f811f 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -41,8 +41,6 @@ #include #include #include -#include -#include #include #include #include @@ -55,36 +53,14 @@ #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 "im-sandbox.h" +#include "process_image.h" /* Write statistics at APPROXIMATELY this interval */ #define STATS_EVERY_N_SECONDS (5) -/* Information about the indexing process for one pattern */ -struct pattern_args -{ - /* "Input" */ - char *filename; - - /* "Output" */ - int n_crystals; -}; - - struct sb_reader { pthread_mutex_t lock; @@ -194,190 +170,6 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead, } -static void process_image(const struct index_args *iargs, - struct pattern_args *pargs, Stream *st, - int cookie) -{ - float *data_for_measurement; - size_t data_size; - int check; - struct hdfile *hdfile; - struct image image; - int i; - char filename[1024]; - - /* Prefix to jump out of temporary folder */ - if ( pargs->filename[0] != '/' ) { - snprintf(filename, 1023, "../../%s", pargs->filename); - } else { - snprintf(filename, 1023, "%s", pargs->filename); - } - - image.features = NULL; - image.data = NULL; - image.flags = NULL; - image.copyme = iargs->copyme; - image.id = cookie; - image.filename = pargs->filename; /* Relative to top level */ - image.beam = iargs->beam; - image.det = iargs->det; - image.crystals = NULL; - image.n_crystals = 0; - - hdfile = hdfile_open(filename); /* Relative to temporary folder */ - if ( hdfile == NULL ) return; - - if ( iargs->element != NULL ) { - - int r; - r = hdfile_set_image(hdfile, iargs->element); - if ( r ) { - ERROR("Couldn't select path '%s'\n", iargs->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; - } - - } - - 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)) - { - 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); - return; - } - - fill_in_values(image.det, hdfile); - fill_in_beam_parameters(image.beam, hdfile); - - image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy)); - - if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) { - /* Error message covers a silly value in the beam file or in - * the HDF5 file. */ - ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n", - image.lambda, image.beam->photon_energy, image.filename); - hdfile_close(hdfile); - return; - } - - if ( iargs->cmfilter ) filter_cm(&image); - - /* Take snapshot of image after CM subtraction but before - * the aggressive noise filter. */ - data_size = image.width * image.height * sizeof(float); - data_for_measurement = malloc(data_size); - - if ( iargs->noisefilter ) { - filter_noise(&image, data_for_measurement); - } else { - memcpy(data_for_measurement, image.data, data_size); - } - - 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"); - } - if ( !iargs->no_revalidate ) { - validate_peaks(&image, iargs->min_int_snr, - iargs->ir_inn, iargs->ir_mid, - iargs->ir_out, iargs->use_saturated); - } - break; - - case PEAK_ZAEF: - search_peaks(&image, iargs->threshold, - iargs->min_gradient, iargs->min_snr, - iargs->ir_inn, iargs->ir_mid, iargs->ir_out, - iargs->use_saturated); - 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; - - /* Index the pattern */ - index_pattern(&image, iargs->indm, iargs->ipriv); - - pargs->n_crystals = image.n_crystals; - - /* Default beam parameters */ - image.div = image.beam->divergence; - image.bw = image.beam->bandwidth; - - /* Integrate each crystal's diffraction spots */ - for ( i=0; iprofile_radius); - - if ( iargs->integrate_found ) { - reflections = select_intersections(&image, - image.crystals[i]); - } else { - reflections = find_intersections(&image, - image.crystals[i]); - } - - crystal_set_reflections(image.crystals[i], reflections); - - } - - /* Integrate all the crystals at once - need all the crystals so that - * overlaps can be detected. */ - integrate_reflections(&image, iargs->closer, - iargs->bgsub, - iargs->min_int_snr, - iargs->ir_inn, - iargs->ir_mid, - iargs->ir_out, - iargs->integrate_saturated, - iargs->res_cutoff); - - write_chunk(st, &image, hdfile, - iargs->stream_peaks, iargs->stream_refls); - - for ( i=0; i + * + * 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 . + * + */ + +#ifdef HAVE_CONFIG_H +#include +#endif + +#include +#include +#include + +#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 "process_image.h" + + +void process_image(const struct index_args *iargs, struct pattern_args *pargs, + Stream *st, int cookie) +{ + float *data_for_measurement; + size_t data_size; + int check; + struct hdfile *hdfile; + struct image image; + int i; + char filename[1024]; + + /* Prefix to jump out of temporary folder */ + if ( pargs->filename[0] != '/' ) { + snprintf(filename, 1023, "../../%s", pargs->filename); + } else { + snprintf(filename, 1023, "%s", pargs->filename); + } + + image.features = NULL; + image.data = NULL; + image.flags = NULL; + image.copyme = iargs->copyme; + image.id = cookie; + image.filename = pargs->filename; /* Relative to top level */ + image.beam = iargs->beam; + image.det = iargs->det; + image.crystals = NULL; + image.n_crystals = 0; + + hdfile = hdfile_open(filename); /* Relative to temporary folder */ + if ( hdfile == NULL ) return; + + if ( iargs->element != NULL ) { + + int r; + r = hdfile_set_image(hdfile, iargs->element); + if ( r ) { + ERROR("Couldn't select path '%s'\n", iargs->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; + } + + } + + 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)) + { + 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); + return; + } + + fill_in_values(image.det, hdfile); + fill_in_beam_parameters(image.beam, hdfile); + + image.lambda = ph_en_to_lambda(eV_to_J(image.beam->photon_energy)); + + if ( (image.beam->photon_energy < 0.0) || (image.lambda > 1000) ) { + /* Error message covers a silly value in the beam file or in + * the HDF5 file. */ + ERROR("Nonsensical wavelength (%e m or %e eV) value for %s.\n", + image.lambda, image.beam->photon_energy, image.filename); + hdfile_close(hdfile); + return; + } + + if ( iargs->cmfilter ) filter_cm(&image); + + /* Take snapshot of image after CM subtraction but before + * the aggressive noise filter. */ + data_size = image.width * image.height * sizeof(float); + data_for_measurement = malloc(data_size); + + if ( iargs->noisefilter ) { + filter_noise(&image, data_for_measurement); + } else { + memcpy(data_for_measurement, image.data, data_size); + } + + 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"); + } + if ( !iargs->no_revalidate ) { + validate_peaks(&image, iargs->min_int_snr, + iargs->ir_inn, iargs->ir_mid, + iargs->ir_out, iargs->use_saturated); + } + break; + + case PEAK_ZAEF: + search_peaks(&image, iargs->threshold, + iargs->min_gradient, iargs->min_snr, + iargs->ir_inn, iargs->ir_mid, iargs->ir_out, + iargs->use_saturated); + 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; + + /* Index the pattern */ + index_pattern(&image, iargs->indm, iargs->ipriv); + + pargs->n_crystals = image.n_crystals; + + /* Default beam parameters */ + image.div = image.beam->divergence; + image.bw = image.beam->bandwidth; + + /* Integrate each crystal's diffraction spots */ + for ( i=0; iprofile_radius); + + if ( iargs->integrate_found ) { + reflections = select_intersections(&image, + image.crystals[i]); + } else { + reflections = find_intersections(&image, + image.crystals[i]); + } + + crystal_set_reflections(image.crystals[i], reflections); + + } + + /* Integrate all the crystals at once - need all the crystals so that + * overlaps can be detected. */ + integrate_reflections(&image, iargs->closer, + iargs->bgsub, + iargs->min_int_snr, + iargs->ir_inn, + iargs->ir_mid, + iargs->ir_out, + iargs->integrate_saturated, + iargs->res_cutoff); + + write_chunk(st, &image, hdfile, + iargs->stream_peaks, iargs->stream_refls); + + for ( i=0; i + * + * 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 . + * + */ + +#ifndef PROCESS_IMAGE_H +#define PROCESS_IMAGE_H + +#ifdef HAVE_CONFIG_H +#include +#endif + + +enum { + PEAK_ZAEF, + PEAK_HDF5, +}; + + +/* Information about the indexing process which is common to all patterns */ +struct index_args +{ + UnitCell *cell; + int cmfilter; + int noisefilter; + int satcorr; + int closer; + int bgsub; + float threshold; + float min_gradient; + float min_snr; + double min_int_snr; + struct detector *det; + IndexingMethod *indm; + IndexingPrivate **ipriv; + int peaks; /* Peak detection method */ + float tols[4]; + struct beam_params *beam; + char *element; + char *hdf5_peak_path; + float ir_inn; + float ir_mid; + float ir_out; + struct copy_hdf5_field *copyme; + int integrate_saturated; + int use_saturated; + int no_revalidate; + int integrate_found; + int stream_peaks; + int stream_refls; + int res_cutoff; +}; + + +/* Information about the indexing process for one pattern */ +struct pattern_args +{ + /* "Input" */ + char *filename; + + /* "Output" */ + int n_crystals; +}; + + +extern void process_image(const struct index_args *iargs, + struct pattern_args *pargs, Stream *st, + int cookie); + + +#endif /* PROCESS_IMAGEs_H */ -- cgit v1.2.3 From fe8352bf60638f0da55669bd9e55a75ad7e300e4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 17:39:26 +0100 Subject: Keep Valgrind quiet --- src/im-sandbox.c | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 290f811f..1ba6b2ce 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -127,6 +127,7 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead, { char *line; char *filename; + size_t len; do { @@ -160,7 +161,13 @@ static char *get_pattern(FILE *fh, char **use_this_one_instead, line = tmp; } - filename = malloc(strlen(prefix)+strlen(line)+1); + len = strlen(prefix)+strlen(line)+1; + + /* Round the length of the buffer, too keep Valgrind quiet when it gets + * given to write() a bit later on */ + len += 4 - (len % 4); + + filename = malloc(len); snprintf(filename, 1023, "%s%s", prefix, line); -- cgit v1.2.3 From 2c439230bad3907f5fc93b3e69ade0d5a5acc4b4 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 18:12:20 +0100 Subject: Avoid the use of stdio for writing to final stream Something really weird goes on when using stdio --- src/im-sandbox.c | 35 +++++++++++++++++++++-------------- src/im-sandbox.h | 2 +- src/indexamajig.c | 11 +++++++---- 3 files changed, 29 insertions(+), 19 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 1ba6b2ce..787e43f9 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -74,8 +74,8 @@ struct sb_reader FILE **fhs; int *fds; - /* Final output file handle */ - FILE *ofh; + /* Final output fd */ + int ofd; }; @@ -277,8 +277,17 @@ static time_t get_monotonic_seconds() #endif +size_t vol = 0; -static int pump_chunk(FILE *fh, FILE *ofh) + +static ssize_t lwrite(int fd, const char *a) +{ + size_t l = strlen(a); + return write(fd, a, l); +} + + +static int pump_chunk(FILE *fh, int ofd) { int chunk_started = 0; @@ -294,25 +303,23 @@ static int pump_chunk(FILE *fh, FILE *ofh) /* Whoops, connection lost */ if ( chunk_started ) { ERROR("EOF during chunk!\n"); - fprintf(ofh, "Chunk is unfinished!\n"); - fprintf(ofh, CHUNK_END_MARKER"\n"); + lwrite(ofd, "Unfinished chunk!\n"); + lwrite(ofd, CHUNK_END_MARKER"\n"); } /* else normal end of output */ return 1; - } else { - ERROR("fgets() failed: %s\n", strerror(errno)); } - break; + + ERROR("fgets() failed: %s\n", strerror(errno)); + return 1; } if ( strcmp(line, "FLUSH\n") == 0 ) break; - - fprintf(ofh, "%s", line); + lwrite(ofd, line); if ( strcmp(line, CHUNK_END_MARKER"\n") == 0 ) break; if ( strcmp(line, CHUNK_START_MARKER"\n") == 0 ) break; - } while ( 1 ); return 0; } @@ -427,7 +434,7 @@ static void *run_reader(void *rdv) /* If the chunk cannot be read, assume the connection * is broken and that the process will die soon. */ - if ( pump_chunk(rd->fhs[i], rd->ofh) ) { + if ( pump_chunk(rd->fhs[i], rd->ofd) ) { /* remove_pipe() assumes that the caller is * holding rd->lock ! */ remove_pipe(rd, i); @@ -632,7 +639,7 @@ static void handle_zombie(struct sandbox *sb) void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, char *use_this_one_instead, - FILE *ofh, int argc, char *argv[]) + int ofd, int argc, char *argv[]) { int i; int allDone; @@ -669,7 +676,7 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, sb->reader->fds = NULL; sb->reader->fhs = NULL; - sb->reader->ofh = ofh; + sb->reader->ofd = ofd; sb->stream_pipe_write = calloc(n_proc, sizeof(int)); if ( sb->stream_pipe_write == NULL ) { diff --git a/src/im-sandbox.h b/src/im-sandbox.h index 7df41738..9248b226 100644 --- a/src/im-sandbox.h +++ b/src/im-sandbox.h @@ -38,5 +38,5 @@ extern void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, int config_basename, FILE *fh, - char *use_this_one_instead, FILE *stream, + char *use_this_one_instead, int streamfd, int argc, char *argv[]); diff --git a/src/indexamajig.c b/src/indexamajig.c index b88469b8..41662d5f 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -44,6 +44,9 @@ #include #include #include +#include +#include +#include #ifdef HAVE_CLOCK_GETTIME #include @@ -158,7 +161,7 @@ int main(int argc, char *argv[]) char *filename = NULL; char *outfile = NULL; FILE *fh; - FILE *ofh; + int ofd; char *rval = NULL; int config_checkprefix = 1; int config_basename = 0; @@ -510,8 +513,8 @@ int main(int argc, char *argv[]) iargs.cell = NULL; } - ofh = fopen(outfile, "w"); - if ( ofh == NULL ) { + ofd = open(outfile, O_CREAT | O_TRUNC | O_WRONLY, S_IRUSR | S_IWUSR); + if ( ofd == -1 ) { ERROR("Failed to open stream '%s'\n", outfile); return 1; } @@ -553,7 +556,7 @@ int main(int argc, char *argv[]) iargs.ipriv = ipriv; create_sandbox(&iargs, n_proc, prefix, config_basename, fh, - use_this_one_instead, ofh, argc, argv); + use_this_one_instead, ofd, argc, argv); free(prefix); -- cgit v1.2.3 From 14dbe0b17405956fd1905044da74c78f09af81b2 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 18:12:58 +0100 Subject: Avoid fclose(NULL) --- src/im-sandbox.c | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index 787e43f9..be96d516 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -543,7 +543,9 @@ static void start_worker_process(struct sandbox *sb, int slot, } for ( j=0; jn_proc; j++ ) { if ( (j != slot) && (sb->running[j]) ) { - fclose(sb->result_fhs[j]); + if ( sb->result_fhs[j] != NULL ) { + fclose(sb->result_fhs[j]); + } close(sb->filename_pipes[j]); } } -- cgit v1.2.3 From 1094e0eeccd8fd68f7a58baea8e03743f89b5edb Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 18:13:11 +0100 Subject: If no progress, still check if everything's finished or not --- src/im-sandbox.c | 2 -- 1 file changed, 2 deletions(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index be96d516..f969e42f 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -785,8 +785,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, break; } - if ( r == 0 ) continue; /* No progress this time. Try again */ - if ( FD_ISSET(signal_pipe[0], &fds) ) { char d; -- cgit v1.2.3 From cbd4ddc97f705b372f3b69f78ba2216d9cdc2d08 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 22 Feb 2013 18:13:33 +0100 Subject: Fussiness --- src/im-sandbox.c | 1 - 1 file changed, 1 deletion(-) (limited to 'src') diff --git a/src/im-sandbox.c b/src/im-sandbox.c index f969e42f..03c553b7 100644 --- a/src/im-sandbox.c +++ b/src/im-sandbox.c @@ -885,7 +885,6 @@ void create_sandbox(struct index_args *iargs, int n_proc, char *prefix, for ( i=0; irunning[i] ) allDone = 0; } - unlock_sandbox(sb); } -- cgit v1.2.3 From ad94c430447aa969397db2356571c62e9bfe4f29 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 27 Feb 2013 14:57:27 +0100 Subject: partial_sim: Initialise image.indexed_by --- src/partial_sim.c | 1 + 1 file changed, 1 insertion(+) (limited to 'src') diff --git a/src/partial_sim.c b/src/partial_sim.c index e48c51fa..a4a2f875 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -530,6 +530,7 @@ int main(int argc, char *argv[]) image.copyme = NULL; image.crystals = NULL; image.n_crystals = 0; + image.indexed_by = INDEXING_NONE; if ( random_intensities ) { full = reflist_new(); -- cgit v1.2.3 From 13659bf4eff299756603920aefef91f227b28186 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 1 Mar 2013 17:02:40 +0100 Subject: indexamajig: Fix default behaviour for --integrate-saturated --- src/indexamajig.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/indexamajig.c b/src/indexamajig.c index 41662d5f..9b868b9f 100644 --- a/src/indexamajig.c +++ b/src/indexamajig.c @@ -204,7 +204,7 @@ int main(int argc, char *argv[]) iargs.ir_mid = 5.0; iargs.ir_out = 7.0; iargs.use_saturated = 0; - iargs.integrate_saturated = 1; + iargs.integrate_saturated = 0; iargs.no_revalidate = 0; iargs.integrate_found = 0; iargs.stream_peaks = 1; -- cgit v1.2.3 From f5f19825380b91edf017af4f7353faf95f97508b Mon Sep 17 00:00:00 2001 From: Thomas White Date: Mon, 4 Mar 2013 10:11:59 +0100 Subject: partialator: Fix various memory bugs --- src/partialator.c | 38 +++++++++++++++++++++++++++++--------- src/scaling-report.c | 7 +++++++ 2 files changed, 36 insertions(+), 9 deletions(-) (limited to 'src') diff --git a/src/partialator.c b/src/partialator.c index 8b51f826..5613c14b 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -410,7 +410,7 @@ int main(int argc, char *argv[]) ERROR("Failed to open input stream '%s'\n", infile); return 1; } - free(infile); + /* Don't free "infile", because it's needed for the scaling report */ /* Sanitise output filename */ if ( outfile == NULL ) { @@ -453,7 +453,6 @@ int main(int argc, char *argv[]) gsl_set_error_handler_off(); /* Fill in what we know about the images so far */ - nobs = 0; n_images = 0; n_crystals = 0; images = NULL; @@ -495,6 +494,7 @@ int main(int argc, char *argv[]) Crystal *cr; Crystal **crystals_new; + RefList *cr_refl; crystals_new = realloc(crystals, (n_crystals+1)*sizeof(Crystal *)); @@ -506,20 +506,20 @@ int main(int argc, char *argv[]) crystals = crystals_new; crystals[n_crystals] = cur->crystals[i]; cr = crystals[n_crystals]; - crystal_set_image(cr, cur); + + /* Image pointer will change due to later reallocs */ + crystal_set_image(cr, NULL); /* Fill in initial estimates of stuff */ crystal_set_osf(cr, 1.0); crystal_set_profile_radius(cr, beam->profile_radius); - crystal_set_user_flag(cr, 0); + crystal_set_user_flag(cr, n_images); /* This is the raw list of reflections */ - as = asymmetric_indices(crystal_get_reflections(cr), - sym); + cr_refl = crystal_get_reflections(cr); + as = asymmetric_indices(cr_refl, sym); crystal_set_reflections(cr, as); - update_partialities(cr); - - nobs += select_scalable_reflections(as, reference); + reflist_free(cr_refl); n_crystals++; @@ -531,6 +531,26 @@ int main(int argc, char *argv[]) close_stream(st); + /* Fill in image pointers */ + nobs = 0; + for ( i=0; i Date: Mon, 4 Mar 2013 10:25:40 +0100 Subject: Free reflists at end --- src/partialator.c | 1 + 1 file changed, 1 insertion(+) (limited to 'src') diff --git a/src/partialator.c b/src/partialator.c index 5613c14b..e67e87c5 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -612,6 +612,7 @@ int main(int argc, char *argv[]) /* Clean up */ for ( i=0; i Date: Mon, 4 Mar 2013 15:51:30 +0100 Subject: partialator: Fix message --- src/partialator.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/partialator.c b/src/partialator.c index e67e87c5..98f24f58 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -605,7 +605,7 @@ int main(int argc, char *argv[]) for ( i=0; i Date: Tue, 5 Mar 2013 10:24:52 +0100 Subject: Update copyright dates --- src/partialator.c | 6 +++--- 1 file changed, 3 insertions(+), 3 deletions(-) (limited to 'src') diff --git a/src/partialator.c b/src/partialator.c index 98f24f58..3c37977b 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -3,11 +3,11 @@ * * Scaling and post refinement for coherent nanocrystallography * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * * This file is part of CrystFEL. * -- cgit v1.2.3 From ba8815d66d9d5a0b39042a756fa04e8d0eda73d6 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 5 Mar 2013 10:28:15 +0100 Subject: partialator: Fix user_flag --- src/partialator.c | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) (limited to 'src') diff --git a/src/partialator.c b/src/partialator.c index 3c37977b..a4af3c18 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -513,7 +513,7 @@ int main(int argc, char *argv[]) /* Fill in initial estimates of stuff */ crystal_set_osf(cr, 1.0); crystal_set_profile_radius(cr, beam->profile_radius); - crystal_set_user_flag(cr, n_images); + crystal_set_user_flag(cr, 0); /* This is the raw list of reflections */ cr_refl = crystal_get_reflections(cr); -- cgit v1.2.3 From 1e03ed982741fdc576ec5a915da120450df20499 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Mar 2013 15:42:23 +0100 Subject: process_hkl: Restore --start-after and --stop-after --- src/process_hkl.c | 45 ++++++++++++++++++++++++++++++++++----------- 1 file changed, 34 insertions(+), 11 deletions(-) (limited to 'src') diff --git a/src/process_hkl.c b/src/process_hkl.c index a8abc2e8..8705ce0f 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -64,8 +64,8 @@ static void show_help(const char *s) " Default: processed.hkl).\n" " -y, --symmetry= Merge according to point group .\n" "\n" -" --start-after= Skip n patterns at the start of the stream.\n" -" --stop-after= Stop after processing n patterns.\n" +" --start-after= Skip crystals at the start of the stream.\n" +" --stop-after= Stop after merging crystals.\n" " -g, --histogram= Calculate the histogram of measurements for this\n" " reflection.\n" " -z, --hist-parameters Set the range for the histogram and the number of\n" @@ -282,7 +282,8 @@ static int merge_all(Stream *st, RefList *model, RefList *reference, const SymOpList *sym, double *hist_vals, signed int hist_h, signed int hist_k, signed int hist_l, - int *hist_i, int config_nopolar, int min_measurements) + int *hist_i, int config_nopolar, int min_measurements, + int start_after, int stop_after) { int rval; int n_images = 0; @@ -290,6 +291,7 @@ static int merge_all(Stream *st, RefList *model, RefList *reference, int n_crystals_used = 0; Reflection *refl; RefListIterator *iter; + int n_crystals_seen = 0; do { @@ -309,6 +311,9 @@ static int merge_all(Stream *st, RefList *model, RefList *reference, int r; Crystal *cr = image.crystals[i]; + n_crystals_seen++; + if ( n_crystals_seen <= start_after ) continue; + n_crystals++; r = merge_crystal(model, &image, cr, reference, sym, hist_vals, hist_h, hist_k, hist_l, @@ -318,12 +323,16 @@ static int merge_all(Stream *st, RefList *model, RefList *reference, crystal_free(cr); + if ( n_crystals_used == stop_after ) break; + } free(image.filename); image_feature_list_free(image.features); - display_progress(n_images, n_crystals, n_crystals_used); + display_progress(n_images, n_crystals_seen, n_crystals_used); + + if ( (stop_after>0) && (n_crystals_used == stop_after) ) break; } while ( rval == 0 ); @@ -377,6 +386,8 @@ int main(int argc, char *argv[]) char *rval; int min_measurements = 2; int r; + int start_after = 0; + int stop_after = 0; /* Long options */ const struct option longopts[] = { @@ -385,8 +396,8 @@ int main(int argc, char *argv[]) {"output", 1, NULL, 'o'}, {"max-only", 0, &config_maxonly, 1}, {"output-every", 1, NULL, 'e'}, - {"stop-after", 1, NULL, 's'}, - {"start-after", 1, NULL, 'f'}, + {"start-after", 1, NULL, 's'}, + {"stop-after", 1, NULL, 'f'}, {"sum", 0, &config_sum, 1}, {"scale", 0, &config_scale, 1}, {"no-polarisation", 0, &config_nopolar, 1}, @@ -399,7 +410,7 @@ int main(int argc, char *argv[]) }; /* Short options */ - while ((c = getopt_long(argc, argv, "hi:e:o:y:g:f:b:z:", + while ((c = getopt_long(argc, argv, "hi:e:o:y:g:s:f:z:", longopts, NULL)) != -1) { switch (c) { @@ -417,11 +428,21 @@ int main(int argc, char *argv[]) break; case 's' : - ERROR("The option '--stop-after' no longer works.\n"); + errno = 0; + start_after = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value for --start-after.\n"); + return 1; + } break; case 'f' : - ERROR("The option '--start-after' no longer works.\n"); + errno = 0; + stop_after = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value for --stop-after.\n"); + return 1; + } break; case 'y' : @@ -525,7 +546,8 @@ int main(int argc, char *argv[]) hist_i = 0; r = merge_all(st, model, NULL, sym, hist_vals, hist_h, hist_k, hist_l, - &hist_i, config_nopolar, min_measurements); + &hist_i, config_nopolar, min_measurements, + start_after, stop_after); fprintf(stderr, "\n"); if ( r ) { ERROR("Error while reading stream.\n"); @@ -554,7 +576,8 @@ int main(int argc, char *argv[]) r = merge_all(st, model, reference, sym, hist_vals, hist_h, hist_k, hist_l, &hist_i, - config_nopolar, min_measurements); + config_nopolar, min_measurements, + start_after, stop_after); fprintf(stderr, "\n"); if ( r ) { ERROR("Error while reading stream.\n"); -- cgit v1.2.3