diff options
| -rw-r--r-- | src/hrs-scaling.c | 14 | ||||
| -rw-r--r-- | src/hrs-scaling.h | 10 | ||||
| -rw-r--r-- | src/partialator.c | 88 |
3 files changed, 93 insertions, 19 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 81c6536b..1c0281a4 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -352,8 +352,8 @@ static void finalise_merge_job(void *vqargs, void *vwargs) } -static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads, - PartialityModel pmodel) +RefList *lsq_intensities(Crystal **crystals, int n, int n_threads, + PartialityModel pmodel) { RefList *full; struct merge_queue_args qargs; @@ -448,9 +448,8 @@ static int test_convergence(double *old_osfs, int n, Crystal **crystals) /* Scale the stack of images */ -RefList *scale_intensities(Crystal **crystals, int n, - int n_threads, int noscale, PartialityModel pmodel, - int min_redundancy) +RefList *scale_intensities(Crystal **crystals, int n, int n_threads, + PartialityModel pmodel, int min_redundancy) { int i; RefList *full = NULL; @@ -463,11 +462,6 @@ RefList *scale_intensities(Crystal **crystals, int n, crystal_set_Bfac(crystals[i], 0.0); } - if ( noscale ) { - full = lsq_intensities(crystals, n, n_threads, pmodel); - return full; - } - /* Create an initial list to refine against */ full = lsq_intensities(crystals, n, n_threads, pmodel); diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h index 16368b79..c41cb323 100644 --- a/src/hrs-scaling.h +++ b/src/hrs-scaling.h @@ -3,11 +3,11 @@ * * Intensity scaling using generalised HRS target function * - * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2013 Thomas White <taw@physics.org> + * 2010-2015 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -40,8 +40,10 @@ #include "geometry.h" extern RefList *scale_intensities(Crystal **crystals, int n, int n_threads, - int noscale, PartialityModel pmodel, - int min_redundancy); + PartialityModel pmodel, int min_redundancy); +extern RefList *lsq_intensities(Crystal **crystals, int n, int n_threads, + PartialityModel pmodel); + #endif /* HRS_SCALING_H */ diff --git a/src/partialator.c b/src/partialator.c index a756ae3d..4b1d5b18 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -244,6 +244,45 @@ static RefList *apply_max_adu(RefList *list, double max_adu) } +static void skip_to_end(FILE *fh) +{ + int c; + do { + c = fgetc(fh); + } while ( (c != '\n') && (c != EOF) ); +} + + +static int set_initial_params(Crystal *cr, FILE *fh) +{ + if ( fh != NULL ) { + + int err; + int n; + float osf, B; + + err = fscanf(fh, "%i %f %f", &n, &osf, &B); + if ( err != 3 ) { + ERROR("Failed to read parameters.\n"); + return 1; + } + + crystal_set_osf(cr, osf); + crystal_set_Bfac(cr, B*1e-20); + + skip_to_end(fh); + + } else { + + crystal_set_osf(cr, 1.0); + crystal_set_Bfac(cr, 0.0); + + } + + return 0; +} + + static void show_duds(Crystal **crystals, int n_crystals) { int j; @@ -326,6 +365,8 @@ int main(int argc, char *argv[]) struct srdata srdata; int polarisation = 1; double max_adu = +INFINITY; + char *sparams_fn = NULL; + FILE *sparams_fh; /* Long options */ const struct option longopts[] = { @@ -341,6 +382,7 @@ int main(int argc, char *argv[]) {"min-measurements", 1, NULL, 2}, {"max-adu", 1, NULL, 3}, + {"start-params", 1, NULL, 4}, {"no-scale", 0, &noscale, 1}, {"no-polarisation", 0, &polarisation, 0}, @@ -415,6 +457,10 @@ int main(int argc, char *argv[]) } break; + case 4 : + sparams_fn = strdup(optarg); + break; + case 0 : break; @@ -473,6 +519,20 @@ int main(int argc, char *argv[]) n_crystals = 0; images = NULL; crystals = NULL; + if ( sparams_fn != NULL ) { + char line[1024]; + sparams_fh = fopen(sparams_fn, "r"); + if ( sparams_fh == NULL ) { + ERROR("Failed to open '%s'\n", sparams_fn); + return 1; + } + fgets(line, 1024, sparams_fh); + STATUS("Reading initial scaling factors (G,B) from '%s'\n", + sparams_fn); + free(sparams_fn); + } else { + sparams_fh = NULL; + } do { @@ -540,6 +600,11 @@ int main(int argc, char *argv[]) crystal_set_user_flag(cr, 0); reflist_free(cr_refl); + if ( set_initial_params(cr, sparams_fh) ) { + ERROR("Failed to set initial parameters\n"); + return 1; + } + n_crystals++; } @@ -551,6 +616,7 @@ int main(int argc, char *argv[]) } while ( 1 ); display_progress(n_images, n_crystals); fprintf(stderr, "\n"); + if ( sparams_fh != NULL ) fclose(sparams_fh); close_stream(st); @@ -581,9 +647,14 @@ int main(int argc, char *argv[]) /* Make initial estimates */ STATUS("Performing initial scaling.\n"); - if ( noscale ) STATUS("Scale factors fixed at 1.\n"); - full = scale_intensities(crystals, n_crystals, - nthreads, noscale, pmodel, min_measurements); + if ( noscale ) { + STATUS("Skipping scaling step (--no-scale).\n"); + full = lsq_intensities(crystals, n_crystals, nthreads, pmodel); + } else { + full = scale_intensities(crystals, n_crystals, nthreads, pmodel, + min_measurements); + } + check_rejection(crystals, n_crystals); srdata.crystals = crystals; @@ -614,8 +685,15 @@ int main(int argc, char *argv[]) /* Re-estimate all the full intensities */ reflist_free(full); - full = scale_intensities(crystals, n_crystals, nthreads, - noscale, pmodel, min_measurements); + if ( noscale ) { + STATUS("Skipping scaling step (--no-scale).\n"); + full = lsq_intensities(crystals, n_crystals, nthreads, + pmodel); + } else { + full = scale_intensities(crystals, n_crystals, nthreads, + pmodel, min_measurements); + } + check_rejection(crystals, n_crystals); srdata.full = full; |