diff options
author | Thomas White <taw@physics.org> | 2018-03-29 16:00:50 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2018-03-29 16:00:50 +0200 |
commit | 8d9f9cc570a69780d3e3b2c271f03d055066a9be (patch) | |
tree | 44d7a71e0cbcab2c1c327d761da150bca9c1807a /src | |
parent | 77cf2edd09bb01ae331935f467064c751f6e338e (diff) | |
parent | 402b3870e3e1bceb974ad5a402c5e5e898f4c87e (diff) |
Merge branch 'tom/ginn-partials'
Diffstat (limited to 'src')
-rw-r--r-- | src/merge.c | 13 | ||||
-rw-r--r-- | src/merge.h | 3 | ||||
-rw-r--r-- | src/partial_sim.c | 5 | ||||
-rw-r--r-- | src/partialator.c | 448 | ||||
-rw-r--r-- | src/post-refinement.c | 1178 | ||||
-rw-r--r-- | src/post-refinement.h | 21 | ||||
-rw-r--r-- | src/process_image.c | 4 | ||||
-rw-r--r-- | src/rejection.c | 2 | ||||
-rw-r--r-- | src/scaling.c | 215 | ||||
-rw-r--r-- | src/scaling.h | 13 |
10 files changed, 1269 insertions, 633 deletions
diff --git a/src/merge.c b/src/merge.c index 17e57c90..e60aca15 100644 --- a/src/merge.c +++ b/src/merge.c @@ -50,7 +50,7 @@ /* Minimum partiality of a reflection for it to be merged */ -#define MIN_PART_MERGE (0.05) +#define MIN_PART_MERGE (0.3) struct merge_queue_args @@ -59,7 +59,6 @@ struct merge_queue_args pthread_rwlock_t full_lock; Crystal **crystals; int n_started; - PartialityModel pmodel; double push_res; int use_weak; long long int n_reflections; @@ -159,7 +158,6 @@ static void run_merge_job(void *vwargs, int cookie) signed int h, k, l; double mean, sumweight, M2, temp, delta, R; double corr, res, w; - //double esd; if ( get_partiality(refl) < MIN_PART_MERGE ) continue; @@ -189,7 +187,7 @@ static void run_merge_job(void *vwargs, int cookie) } /* Total (multiplicative) correction factor */ - corr = exp(-G) * exp(B*res*res) * get_lorentz(refl) + corr = G * exp(B*res*res) * get_lorentz(refl) / get_partiality(refl); if ( isnan(corr) ) { ERROR("NaN corr:\n"); @@ -200,8 +198,8 @@ static void run_merge_job(void *vwargs, int cookie) continue; } - //esd = get_esd_intensity(refl) * corr; - w = 1.0; + /* Reflections count less the more they have to be scaled up */ + w = 1.0 / corr; /* Running mean and variance calculation */ temp = w + sumweight; @@ -229,7 +227,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs) RefList *merge_intensities(Crystal **crystals, int n, int n_threads, - PartialityModel pmodel, int min_meas, + int min_meas, double push_res, int use_weak) { RefList *full; @@ -245,7 +243,6 @@ RefList *merge_intensities(Crystal **crystals, int n, int n_threads, qargs.full = full; qargs.n_started = 0; qargs.crystals = crystals; - qargs.pmodel = pmodel; qargs.push_res = push_res; qargs.use_weak = use_weak; qargs.n_reflections = 0; diff --git a/src/merge.h b/src/merge.h index 1911584f..56229e4c 100644 --- a/src/merge.h +++ b/src/merge.h @@ -40,7 +40,6 @@ #include "geometry.h" extern RefList *merge_intensities(Crystal **crystals, int n, int n_threads, - PartialityModel pmodel, int min_meas, - double push_res, int use_weak); + int min_meas, double push_res, int use_weak); #endif /* MERGE */ diff --git a/src/partial_sim.c b/src/partial_sim.c index 4a348742..fc7845c1 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -386,7 +386,7 @@ static void run_job(void *vwargs, int cookie) reflections = predict_to_res(cr, largest_q(&wargs->image)); crystal_set_reflections(cr, reflections); - calculate_partialities(cr, PMODEL_SCSPHERE); + calculate_partialities(cr, PMODEL_XSPHERE); for ( i=0; i<NBINS; i++ ) { wargs->n_ref[i] = 0; @@ -862,7 +862,7 @@ int main(int argc, char *argv[]) STATUS(" Background: none (no image " "output)\n"); } - STATUS(" Partiality model: scsphere (hardcoded)\n"); + STATUS(" Partiality model: xsphere (hardcoded)\n"); STATUS(" Noise standard deviation: %.2f detector units\n", noise_stddev); if ( random_intensities ) { @@ -874,6 +874,7 @@ int main(int argc, char *argv[]) input_file, sym_str); } STATUS(" Max error in cell components: %.2f %%\n", cnoise); + STATUS("Scale factor standard deviation: %.2f\n", osf_stddev); if ( random_intensities ) { full = reflist_new(); diff --git a/src/partialator.c b/src/partialator.c index 428b0d68..21b02266 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -3,11 +3,11 @@ * * Scaling and post refinement for coherent nanocrystallography * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -40,6 +40,7 @@ #include <assert.h> #include <pthread.h> #include <gsl/gsl_errno.h> +#include <sys/stat.h> #include <image.h> #include <utils.h> @@ -190,7 +191,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile, } snprintf(tmp, 1024, "%s1", outfile); split = merge_intensities(crystals1, n_crystals1, nthreads, - pmodel, min_measurements, push_res, 1); + min_measurements, push_res, 1); if ( split == NULL ) { ERROR("Not enough crystals for two way split!\n"); @@ -202,7 +203,7 @@ static void write_split(Crystal **crystals, int n_crystals, const char *outfile, reflist_free(split); snprintf(tmp, 1024, "%s2", outfile); split = merge_intensities(crystals2, n_crystals2, nthreads, - pmodel, min_measurements, push_res, 1); + min_measurements, push_res, 1); STATUS("and %s\n", tmp); write_reflist_2(tmp, split, sym); reflist_free(split); @@ -289,7 +290,7 @@ static void write_custom_split(struct custom_split *csplit, int dsn, STATUS("Writing dataset '%s' to %s (%i crystals)\n", csplit->dataset_names[dsn], tmp, n_crystalsn); split = merge_intensities(crystalsn, n_crystalsn, nthreads, - pmodel, min_measurements, push_res, 1); + min_measurements, push_res, 1); write_reflist_2(tmp, split, sym); reflist_free(split); @@ -325,7 +326,12 @@ static void show_help(const char *s) " -j <n> Run <n> analyses in parallel.\n" " --no-free Disable cross-validation (testing only).\n" " --custom-split List of files for custom dataset splitting.\n" -" --max-rel-B Maximum allowable relative |B| factor.\n"); +" --max-rel-B Maximum allowable relative |B| factor.\n" +" --no-logs Do not write extensive log files.\n" +" -w <pg> Apparent point group for resolving ambiguities.\n" +" --operator=<op> Indexing ambiguity operator for resolving.\n" +" --force-bandwidth=<n> Set all bandwidths to <n> (fraction).\n" +" --force-radius=<n> Set all profile radii to <n> nm^-1.\n"); } @@ -579,7 +585,7 @@ static int set_initial_params(Crystal *cr, FILE *fh) } else { - crystal_set_osf(cr, 0.0); + crystal_set_osf(cr, 1.0); crystal_set_Bfac(cr, 0.0); } @@ -604,13 +610,22 @@ static void select_free_reflections(RefList *list, gsl_rng *rng) static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr, - int fr) + int fr, signed int inum) { Reflection *refl; RefListIterator *iter; double G = crystal_get_osf(cr); double B = crystal_get_Bfac(cr); UnitCell *cell = crystal_get_cell(cr); + char ins[5]; + + if ( inum >= 0 ) { + snprintf(ins, 4, "%i", inum); + } else { + ins[0] = 'F'; + ins[1] = '\0'; + } + for ( refl = first_refl(list, &iter); refl != NULL; @@ -623,6 +638,9 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr, if ( !get_flag(refl) ) continue; /* Not free-flagged */ + /* Strong reflections only */ + if ( get_intensity(refl) < 3.0*get_esd_intensity(refl) ) continue; + get_indices(refl, &h, &k, &l); res = resolution(cell, h, k, l); if ( 2.0*res > crystal_get_resolution_limit(cr) ) continue; @@ -630,46 +648,55 @@ static void write_to_pgraph(FILE *fh, RefList *list, RefList *full, Crystal *cr, match = find_refl(full, h, k, l); if ( match == NULL ) continue; + /* Don't calculate pobs if reference reflection is weak */ + if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue; + /* Calculated partiality */ pcalc = get_partiality(refl); /* Observed partiality */ - corr = exp(-G) * exp(B*res*res) * get_lorentz(refl); + corr = G * exp(B*res*res) * get_lorentz(refl); Ipart = get_intensity(refl) * corr; pobs = Ipart / get_intensity(match); - fprintf(fh, "%5i %4i %4i %4i %8.4f %8.3f %8.3f\n", - fr, h, k, l, 2*res/1e9, pcalc, pobs); + fprintf(fh, "%5i %4i %4i %4i %e %e %8.3f %8.3f %s\n", + fr, h, k, l, 2*res, Ipart, pcalc, pobs, ins); } } static void write_pgraph(RefList *full, Crystal **crystals, int n_crystals, - int iter) + signed int iter, const char *suff) { FILE *fh; char tmp[256]; int i; - snprintf(tmp, 256, "pgraph-iter%i.dat", iter); + snprintf(tmp, 256, "pr-logs/pgraph%s.dat", suff); + + if ( iter == 0 ) { + fh = fopen(tmp, "w"); + } else { + fh = fopen(tmp, "a"); + } - fh = fopen(tmp, "w"); if ( fh == NULL ) { ERROR("Failed to open '%s'\n", tmp); return; } - fprintf(fh, " fr h k l 1/d(nm) pcalc pobs\n"); + if ( iter == 0 ) { + fprintf(fh, " Crystal h k l 1/d(m) Ipart pcalc pobs iteration\n"); + } for ( i=0; i<n_crystals; i++ ) { if ( crystal_get_user_flag(crystals[i]) != 0 ) continue; write_to_pgraph(fh, crystal_get_reflections(crystals[i]), full, - crystals[i], i); + crystals[i], i, iter); } fclose(fh); - } @@ -690,8 +717,8 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full, if ( crystal_get_user_flag(crystals[i]) ) continue; - r = residual(crystals[i], full, 0, NULL, NULL); - free_r = residual(crystals[i], full, 1, NULL, NULL); + r = residual(crystals[i], full, 0, NULL, NULL, 1); + free_r = residual(crystals[i], full, 1, NULL, NULL, 1); log_r = log_residual(crystals[i], full, 0, NULL, NULL); free_log_r = log_residual(crystals[i], full, 1, NULL, NULL); @@ -720,32 +747,77 @@ static void show_all_residuals(Crystal **crystals, int n_crystals, } -static void dump_parameters(const char *filename, Crystal **crystals, - int n_crystals) +struct log_qargs { - FILE *fh; - fh = fopen(filename, "w"); - if ( fh == NULL ) { - ERROR("Couldn't open partialator.params!\n"); - } else { - fprintf(fh, " cr OSF relB div" - " flag filename event\n"); - int i; - for ( i=0; i<n_crystals; i++ ) { - struct image *img; - char *evt_str; - img = crystal_get_image(crystals[i]); - evt_str = get_event_string(img->event); - fprintf(fh, "%4i %10.5f %10.2f %8.5e %-25s %s %s\n", - i, crystal_get_osf(crystals[i]), - crystal_get_Bfac(crystals[i])*1e20, - crystal_get_image(crystals[i])->div, - str_prflag(crystal_get_user_flag(crystals[i])), - img->filename, evt_str); - free(evt_str); - } - fclose(fh); - } + int iter; + int next; + Crystal **crystals; + int n_crystals; + RefList *full; +}; + + +struct log_args +{ + Crystal *cr; + RefList *full; + int iter; + int cnum; +}; + + +static void *get_log_task(void *vp) +{ + struct log_qargs *qargs = vp; + struct log_args *task; + + if ( qargs->next >= qargs->n_crystals ) return NULL; + + task = malloc(sizeof(struct log_args)); + if ( task == NULL ) return NULL; + + task->cr = qargs->crystals[qargs->next]; + task->full = qargs->full; + task->iter = qargs->iter; + task->cnum = qargs->next; + + qargs->next += 20; + return task; +} + + +static void write_logs(void *vp, int cookie) +{ + struct log_args *args = vp; + write_specgraph(args->cr, args->full, args->iter, args->cnum); + write_gridscan(args->cr, args->full, args->iter, args->cnum); + write_test_logs(args->cr, args->full, args->iter, args->cnum); +} + + +static void done_log(void *qargs, void *vp) +{ + struct log_args *task = vp; + free(task); +} + + +static void write_logs_parallel(Crystal **crystals, int n_crystals, + RefList *full, int iter, int n_threads) +{ + struct log_qargs qargs; + + + qargs.iter = iter; + qargs.next = 0; + qargs.full = full; + qargs.crystals = crystals; + qargs.n_crystals = n_crystals; + + STATUS("Writing logs...\n"); + run_threads(n_threads, write_logs, get_log_task, done_log, &qargs, + n_crystals/20, 0, 0, 0); + STATUS("Done.\n"); } @@ -756,9 +828,10 @@ int main(int argc, char *argv[]) char *outfile = NULL; char *sym_str = NULL; SymOpList *sym; + SymOpList *amb; + SymOpList *w_sym; int nthreads = 1; int i; - struct image *images; int n_iter = 10; RefList *full; int n_images = 0; @@ -770,7 +843,7 @@ int main(int argc, char *argv[]) Stream *st; Crystal **crystals; char *pmodel_str = NULL; - PartialityModel pmodel = PMODEL_SCSPHERE; + PartialityModel pmodel = PMODEL_XSPHERE; int min_measurements = 2; char *rval; int polarisation = 1; @@ -786,6 +859,13 @@ int main(int argc, char *argv[]) char *csplit_fn = NULL; struct custom_split *csplit = NULL; double max_B = 1e-18; + char *rfile = NULL; + RefList *reference = NULL; + int no_logs = 0; + char *w_sym_str = NULL; + char *operator = NULL; + double force_bandwidth = -1.0; + double force_radius = -1.0; /* Long options */ const struct option longopts[] = { @@ -806,8 +886,12 @@ int main(int argc, char *argv[]) {"push-res", 1, NULL, 5}, {"res-push", 1, NULL, 5}, /* compat */ {"custom-split", 1, NULL, 6}, - {"max-rel-B" , 1, NULL, 7}, - {"max-rel-b" , 1, NULL, 7}, /* compat */ + {"max-rel-B", 1, NULL, 7}, + {"max-rel-b", 1, NULL, 7}, /* compat */ + {"reference", 1, NULL, 8}, /* ssshhh! */ + {"operator", 1, NULL, 9}, + {"force-bandwidth", 1, NULL, 10}, + {"force-radius", 1, NULL, 11}, {"no-scale", 0, &no_scale, 1}, {"no-pr", 0, &no_pr, 1}, @@ -817,6 +901,7 @@ int main(int argc, char *argv[]) {"polarization", 0, &polarisation, 1}, /* compat */ {"no-free", 0, &no_free, 1}, {"output-every-cycle", 0, &output_everycycle, 1}, + {"no-logs", 0, &no_logs, 1}, {0, 0, NULL, 0} }; @@ -828,7 +913,7 @@ int main(int argc, char *argv[]) } /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:", + while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:j:m:w:", longopts, NULL)) != -1) { @@ -887,6 +972,10 @@ int main(int argc, char *argv[]) pmodel_str = strdup(optarg); break; + case 'w' : + w_sym_str = strdup(optarg); + break; + case 2 : errno = 0; min_measurements = strtod(optarg, &rval); @@ -933,6 +1022,33 @@ int main(int argc, char *argv[]) max_B = max_B * 1e-20; break; + case 8 : + rfile = strdup(optarg); + break; + + case 9 : + operator = strdup(optarg); + break; + + case 10 : + errno = 0; + force_bandwidth = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value for --force-bandwidth.\n"); + return 1; + } + break; + + case 11 : + errno = 0; + force_radius = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value for --force-radius.\n"); + return 1; + } + force_radius *= 1e9; + break; + case 0 : break; @@ -972,13 +1088,49 @@ int main(int argc, char *argv[]) free(sym_str); if ( sym == NULL ) return 1; + if ( (w_sym_str != NULL) && (operator != NULL) ) { + ERROR("Specify the apparent symmetry (-w) or the operator, " + "not both.\n"); + return 1; + } + + if ( w_sym_str == NULL ) { + w_sym = NULL; + amb = NULL; + } else { + pointgroup_warning(w_sym_str); + w_sym = get_pointgroup(w_sym_str); + free(w_sym_str); + if ( w_sym == NULL ) return 1; + amb = get_ambiguities(w_sym, sym); + if ( amb == NULL ) { + ERROR("Couldn't find ambiguity operator.\n"); + ERROR("Check that your values for -y and -w are " + "correct.\n"); + return 1; + } + + } + + if ( operator ) { + amb = parse_symmetry_operations(operator); + if ( amb == NULL ) return 1; + set_symmetry_name(amb, "Ambiguity"); + } + + if ( amb != NULL ) { + STATUS("Indexing ambiguity resolution enabled. " + "The ambiguity operation(s) are:\n"); + describe_symmetry(amb); + /* In contrast to ambigator, partialator can deal with multiple + * ambiguities at once */ + } + if ( pmodel_str != NULL ) { if ( strcmp(pmodel_str, "unity") == 0 ) { pmodel = PMODEL_UNITY; - } else if ( strcmp(pmodel_str, "scgaussian") == 0 ) { - pmodel = PMODEL_SCGAUSSIAN; - } else if ( strcmp(pmodel_str, "scsphere") == 0 ) { - pmodel = PMODEL_SCSPHERE; + } else if ( strcmp(pmodel_str, "xsphere") == 0 ) { + pmodel = PMODEL_XSPHERE; } else if ( strcmp(pmodel_str, "random") == 0 ) { pmodel = PMODEL_RANDOM; } else { @@ -989,7 +1141,22 @@ int main(int argc, char *argv[]) if ( (pmodel == PMODEL_UNITY) && !no_pr ) { no_pr = 1; - STATUS("Setting --no-pr because we are not modelling partialities (--model=unity).\n"); + STATUS("Setting --no-pr because we are not modelling " + "partialities (--model=unity).\n"); + } + + if ( !no_logs ) { + int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH); + if ( r ) { + if ( errno == EEXIST ) { + ERROR("A folder called 'pr-logs' exists in the " + "working directory.\n"); + ERROR("Please delete or move it first.\n"); + } else { + ERROR("Failed to create pr-logs folder.\n"); + } + return 1; + } } /* Read the custom split list (if applicable) */ @@ -1002,6 +1169,20 @@ int main(int argc, char *argv[]) free(csplit_fn); } + if ( rfile != NULL ) { + RefList *rread; + rread = read_reflections(rfile); + if ( rread == NULL ) { + ERROR("Failed to read reference reflections\n"); + return 1; + } + reference = asymmetric_indices(rread, sym); + reflist_free(rread); + ERROR("WARNING: Using an external reference.\n"); + ERROR("WARNING: If you publish a structure based on the result," + " expect to have to retract your paper!\n"); + } + gsl_set_error_handler_off(); rng = gsl_rng_alloc(gsl_rng_mt19937); @@ -1009,7 +1190,6 @@ int main(int argc, char *argv[]) n_images = 0; n_crystals = 0; n_crystals_seen = 0; - images = NULL; crystals = NULL; if ( sparams_fn != NULL ) { char line[1024]; @@ -1033,35 +1213,27 @@ int main(int argc, char *argv[]) RefList *as; int i; - struct image *images_new; - struct image *cur; + struct image cur; - images_new = realloc(images, (n_images+1)*sizeof(struct image)); - if ( images_new == NULL ) { - ERROR("Failed to allocate memory for image list.\n"); - return 1; - } - images = images_new; - cur = &images[n_images]; - - cur->div = NAN; - cur->bw = NAN; - cur->det = NULL; - if ( read_chunk_2(st, cur, STREAM_READ_REFLECTIONS - | STREAM_READ_UNITCELL) != 0 ) { + cur.div = NAN; + cur.bw = NAN; + cur.det = NULL; + if ( read_chunk_2(st, &cur, STREAM_READ_REFLECTIONS + | STREAM_READ_UNITCELL) != 0 ) { break; } - if ( isnan(cur->div) || isnan(cur->bw) ) { + if ( isnan(cur.div) || isnan(cur.bw) ) { ERROR("Chunk doesn't contain beam parameters.\n"); return 1; } - for ( i=0; i<cur->n_crystals; i++ ) { + for ( i=0; i<cur.n_crystals; i++ ) { Crystal *cr; Crystal **crystals_new; RefList *cr_refl; + struct image *image; n_crystals_seen++; if ( n_crystals_seen <= start_after ) continue; @@ -1074,11 +1246,19 @@ int main(int argc, char *argv[]) return 1; } crystals = crystals_new; - crystals[n_crystals] = cur->crystals[i]; + crystals[n_crystals] = cur.crystals[i]; cr = crystals[n_crystals]; - /* Image pointer will change due to later reallocs */ - crystal_set_image(cr, NULL); + image = malloc(sizeof(struct image)); + if ( image == NULL ) { + ERROR("Failed to allocatea memory for image.\n"); + return 1; + } + + crystal_set_image(cr, image); + *image = cur; + image->n_crystals = 1; + image->crystals = &crystals[n_crystals]; /* This is the raw list of reflections */ cr_refl = crystal_get_reflections(cr); @@ -1088,7 +1268,7 @@ int main(int argc, char *argv[]) if ( polarisation ) { polarisation_correction(cr_refl, crystal_get_cell(cr), - cur); + image); } if ( !no_free ) select_free_reflections(cr_refl, rng); @@ -1124,65 +1304,76 @@ int main(int argc, char *argv[]) close_stream(st); - /* Fill in image pointers */ - for ( i=0; i<n_images; i++ ) { - int j; - for ( j=0; j<images[i].n_crystals; j++ ) { - - Crystal *cryst; - - cryst = images[i].crystals[j]; - crystal_set_image(cryst, &images[i]); - - /* Now it's safe to do the following */ - update_predictions(cryst); - calculate_partialities(cryst, pmodel); - + STATUS("Initial partiality calculation...\n"); + for ( i=0; i<n_crystals; i++ ) { + Crystal *cr = crystals[i]; + if ( force_bandwidth > 0.0 ) { + crystal_get_image(cr)->bw = force_bandwidth; + } + if ( force_radius > 0.0 ) { + crystal_set_profile_radius(cr, force_radius); } + update_predictions(cr); + calculate_partialities(cr, pmodel); } if (csplit != NULL) check_csplit(crystals, n_crystals, csplit); /* Make a first pass at cutting out crap */ - STATUS("Checking patterns.\n"); + STATUS("Early rejection...\n"); //early_rejection(crystals, n_crystals); - /* Initial rejection, figures of merit etc */ - full = merge_intensities(crystals, n_crystals, nthreads, pmodel, - min_measurements, push_res, 1); + /* Create reference data set if we don't already have one */ + if ( reference == NULL ) { + STATUS("Initial scaling...\n"); + scale_all(crystals, n_crystals, nthreads); + full = merge_intensities(crystals, n_crystals, nthreads, + min_measurements, push_res, 1); + } else { + full = reference; + } + + /* Scale everything to the reference */ + STATUS("Initial scaling to reference...\n"); + if ( !no_scale ) { + scale_all_to_reference(crystals, n_crystals, full, nthreads); + } + + /* Check rejection and write figures of merit */ check_rejection(crystals, n_crystals, full, max_B); show_all_residuals(crystals, n_crystals, full); - write_pgraph(full, crystals, n_crystals, 0); + if ( !no_pr && !no_logs ) { + write_pgraph(full, crystals, n_crystals, 0, ""); + write_logs_parallel(crystals, n_crystals, full, 0, nthreads); + } /* Iterate */ for ( i=0; i<n_iter; i++ ) { STATUS("Scaling and refinement cycle %i of %i\n", i+1, n_iter); - if ( !no_scale ) { - scale_all(crystals, n_crystals, nthreads, pmodel); - reflist_free(full); - full = merge_intensities(crystals, n_crystals, nthreads, - pmodel, min_measurements, - push_res, 1); + if ( !no_pr ) { + refine_all(crystals, n_crystals, full, nthreads, pmodel, + 0, i+1, no_logs, sym, amb); + } else if ( !no_scale ) { + scale_all_to_reference(crystals, n_crystals, full, nthreads); } - if ( !no_pr ) { - refine_all(crystals, n_crystals, full, nthreads, - pmodel); + /* Create new reference if needed */ + if ( reference == NULL ) { reflist_free(full); + scale_all(crystals, n_crystals, nthreads); full = merge_intensities(crystals, n_crystals, nthreads, - pmodel, min_measurements, + min_measurements, push_res, 1); - } + } /* else full still equals reference */ check_rejection(crystals, n_crystals, full, max_B); - reflist_free(full); - full = merge_intensities(crystals, n_crystals, nthreads, - pmodel, min_measurements, - push_res, 1); show_all_residuals(crystals, n_crystals, full); - write_pgraph(full, crystals, n_crystals, i+1); + + if ( !no_logs ) { + write_pgraph(full, crystals, n_crystals, i+1, ""); + } if ( output_everycycle ) { @@ -1209,17 +1400,29 @@ int main(int argc, char *argv[]) } } - /* Dump parameters */ - snprintf(tmp, 1024, "iter%.2d_partialator.params", i+1); - dump_parameters(tmp, crystals, n_crystals); - } } - full = merge_intensities(crystals, n_crystals, nthreads, pmodel, - min_measurements, push_res, 1); + /* Final merge */ + STATUS("Final merge...\n"); + if ( reference == NULL ) { + reflist_free(full); + scale_all(crystals, n_crystals, nthreads); + full = merge_intensities(crystals, n_crystals, nthreads, + min_measurements, + push_res, 1); + } else { + scale_all_to_reference(crystals, n_crystals, reference, nthreads); + full = merge_intensities(crystals, n_crystals, nthreads, + min_measurements, push_res, 1); + } + + /* Write final figures of merit (no rejection any more) */ show_all_residuals(crystals, n_crystals, full); - write_pgraph(full, crystals, n_crystals, n_iter+1); + if ( !no_pr && !no_logs ) { + write_pgraph(full, crystals, n_crystals, -1, ""); + write_logs_parallel(crystals, n_crystals, full, -1, nthreads); + } /* Output results */ STATUS("Writing overall results to %s\n", outfile); @@ -1239,9 +1442,6 @@ int main(int argc, char *argv[]) } } - /* Dump parameters */ - dump_parameters("partialator.params", crystals, n_crystals); - /* Clean up */ gsl_rng_free(rng); for ( i=0; i<n_crystals; i++ ) { @@ -1249,13 +1449,9 @@ int main(int argc, char *argv[]) crystal_free(crystals[i]); } reflist_free(full); - free(sym); + free_symoplist(sym); free(outfile); free(crystals); - for ( i=0; i<n_images; i++ ) { - free(images[i].filename); - } - free(images); free(infile); return 0; diff --git a/src/post-refinement.c b/src/post-refinement.c index bf30d299..79ea187a 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -33,11 +33,8 @@ #include <stdlib.h> #include <assert.h> -#include <gsl/gsl_matrix.h> -#include <gsl/gsl_vector.h> -#include <gsl/gsl_linalg.h> -#include <gsl/gsl_eigen.h> -#include <gsl/gsl_blas.h> +#include <gsl/gsl_multimin.h> +#include <gsl/gsl_fit.h> #include "image.h" #include "post-refinement.h" @@ -46,15 +43,13 @@ #include "geometry.h" #include "cell.h" #include "cell-utils.h" +#include "reflist-utils.h" +#include "scaling.h" -/* Maximum number of iterations of NLSq to do for each image per macrocycle. */ -#define MAX_CYCLES (10) - struct prdata { int refined; - int n_filtered; }; const char *str_prflag(enum prflag flag) @@ -79,625 +74,867 @@ const char *str_prflag(enum prflag flag) case PRFLAG_BIGB : return "B too big"; + case PRFLAG_SCALEBAD : + return "bad scaling"; + default : return "Unknown flag"; } } -/* Returns dp(gauss)/dr at "r" */ -static double gaussian_fraction_gradient(double r, double R) +double residual(Crystal *cr, const RefList *full, int free, + int *pn_used, const char *filename, int complain) { - const double ng = 2.6; - const double sig = R/ng; + Reflection *refl; + RefListIterator *iter; + int n_used = 0; + double num = 0.0; + double den = 0.0; + double G = crystal_get_osf(cr); + double B = crystal_get_Bfac(cr); + UnitCell *cell = crystal_get_cell(cr); - /* If the Ewald sphere isn't within the profile, the gradient is zero */ - if ( r < -R ) return 0.0; - if ( r > +R ) return 0.0; + if ( linear_scale(full, crystal_get_reflections(cr), &G, complain) ) { + return GSL_NAN; + } - return exp(-pow(r/sig, 2.0)/2.0) / (sig*sqrt(2.0*M_PI)); -} + for ( refl = first_refl(crystal_get_reflections(cr), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + double p, w, corr, res; + signed int h, k, l; + Reflection *match; + double I_full; + double int1, int2; + if ( free && !get_flag(refl) ) continue; -/* Returns dp(sph)/dr at "r" */ -static double sphere_fraction_gradient(double r, double pr) -{ - double q, dpdq, dqdr; + get_indices(refl, &h, &k, &l); + res = resolution(cell, h, k, l); + match = find_refl(full, h, k, l); + if ( match == NULL ) continue; + I_full = get_intensity(match); - /* If the Ewald sphere isn't within the profile, the gradient is zero */ - if ( r < -pr ) return 0.0; - if ( r > +pr ) return 0.0; + if ( get_redundancy(match) < 2 ) continue; - q = (r + pr)/(2.0*pr); - dpdq = 6.0*(q - q*q); - dqdr = 1.0 / (2.0*pr); - return dpdq * dqdr; -} + p = get_partiality(refl); + //if ( p < 0.2 ) continue; + corr = G * exp(B*res*res) * get_lorentz(refl); + int1 = get_intensity(refl) * corr; + int2 = p*I_full; + w = p; -/* Returns dp/dr at "r" */ -static double partiality_gradient(double r, double pr, - PartialityModel pmodel, - double rlow, double rhigh) -{ - double A, D; + num += fabs(int1 - int2) * w; + den += fabs(int1 + int2) * w/2.0; - D = rlow - rhigh; + n_used++; - switch ( pmodel ) { + } + + if ( pn_used != NULL ) *pn_used = n_used; + return num/(den*sqrt(2)); +} - default: - case PMODEL_UNITY: - return 0.0; - case PMODEL_SCSPHERE: - A = sphere_fraction_gradient(r, pr)/D; - return 4.0*pr*A/3.0; +static UnitCell *rotate_cell_xy(const UnitCell *cell, double ang1, double ang2) +{ + UnitCell *o; + double asx, asy, asz; + double bsx, bsy, bsz; + double csx, csy, csz; + double xnew, ynew, znew; + + o = cell_new_from_cell(cell); + + cell_get_reciprocal(o, &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + + /* "a" around x */ + xnew = asx; + ynew = asy*cos(ang1) + asz*sin(ang1); + znew = -asy*sin(ang1) + asz*cos(ang1); + asx = xnew; asy = ynew; asz = znew; + + /* "b" around x */ + xnew = bsx; + ynew = bsy*cos(ang1) + bsz*sin(ang1); + znew = -bsy*sin(ang1) + bsz*cos(ang1); + bsx = xnew; bsy = ynew; bsz = znew; + + /* "c" around x */ + xnew = csx; + ynew = csy*cos(ang1) + csz*sin(ang1); + znew = -csy*sin(ang1) + csz*cos(ang1); + csx = xnew; csy = ynew; csz = znew; + + /* "a" around y */ + xnew = asx*cos(ang2) + asz*sin(ang2); + ynew = asy; + znew = -asx*sin(ang2) + asz*cos(ang2); + asx = xnew; asy = ynew; asz = znew; + + /* "b" around y */ + xnew = bsx*cos(ang2) + bsz*sin(ang2); + ynew = bsy; + znew = -bsx*sin(ang2) + bsz*cos(ang2); + bsx = xnew; bsy = ynew; bsz = znew; + + /* "c" around y */ + xnew = csx*cos(ang2) + csz*sin(ang2); + ynew = csy; + znew = -csx*sin(ang2) + csz*cos(ang2); + csx = xnew; csy = ynew; csz = znew; + + cell_set_reciprocal(o, asx, asy, asz, bsx, bsy, bsz, csx, csy, csz); + + return o; +} - case PMODEL_SCGAUSSIAN: - A = gaussian_fraction_gradient(r, pr)/D; - return 4.0*pr*A/3.0; +static const char *get_label(enum gparam p) +{ + switch ( p ) { + case GPARAM_ANG1 : return "First angle (radians)"; + case GPARAM_ANG2 : return "Second angle (radians)"; + case GPARAM_R : return "Profile radius (m^-1)"; + case GPARAM_WAVELENGTH : return "Wavelength (Angstroms)"; + default : return "unknown"; } } -static double sphere_fraction_rgradient(double r, double R) +/* We set all the step sizes to 1, then scale them. + * This way, the size of the simplex stays meaningful and we possibly also + * avoid some roundoff errors */ +static double get_scale(enum gparam p) { - /* If the Ewald sphere isn't within the profile, the gradient is zero */ - if ( r < -R ) return 0.0; - if ( r > +R ) return 0.0; + switch ( p ) { + + case GPARAM_ANG1 : return deg2rad(0.05); + case GPARAM_ANG2 : return deg2rad(0.05); + case GPARAM_R : return 0.0005e9; + case GPARAM_WAVELENGTH : return 0.001e-10; - return -(3.0*r/(4.0*R*R)) * (1.0 - r*r/(R*R)); + default : return 0.0; + + } } -static double gaussian_fraction_rgradient(double r, double R) -{ - const double ng = 2.6; - const double sig = R/ng; +struct rf_priv { + const Crystal *cr; + const RefList *full; + enum gparam rv[32]; + int verbose; + int serial; + gsl_vector *initial; - /* If the Ewald sphere isn't within the profile, the gradient is zero */ - if ( r < -R ) return 0.0; - if ( r > +R ) return 0.0; + /* For freeing later */ + gsl_vector *vals; + gsl_vector *step; - return -(ng*r/(sqrt(2.0*M_PI)*R*R))*exp(-r*r/(2.0*sig*sig)); + /* So that it stays around until the end of minimisation */ + gsl_multimin_function f; +}; + + +static double get_actual_val(const gsl_vector *v, const gsl_vector *initial, + enum gparam *rv, int i) +{ + return gsl_vector_get(v, i) * get_scale(rv[i]) + + gsl_vector_get(initial, i); } -static double volume_fraction_rgradient(double r, double pr, - PartialityModel pmodel) +static void apply_parameters(const gsl_vector *v, const gsl_vector *initial, + enum gparam *rv, Crystal *cr) { - switch ( pmodel ) - { - case PMODEL_UNITY : - return 1.0; + int i; + double ang1, ang2, R, lambda; + UnitCell *cell; - case PMODEL_SCSPHERE : - return sphere_fraction_rgradient(r, pr); + /* Default parameters if not used in refinement */ + ang1 = 0.0; + ang2 = 0.0; + R = crystal_get_profile_radius(cr); + lambda = crystal_get_image(cr)->lambda; - case PMODEL_SCGAUSSIAN : - return gaussian_fraction_rgradient(r, pr); + for ( i=0; i<v->size; i++ ) { - default : - ERROR("No pmodel in volume_fraction_rgradient!\n"); - return 1.0; - } -} + double val; + val = get_actual_val(v, initial, rv, i); -static double volume_fraction(double rlow, double rhigh, double pr, - PartialityModel pmodel) -{ - switch ( pmodel ) - { - case PMODEL_UNITY : - return 1.0; + switch ( rv[i] ) { - case PMODEL_SCSPHERE : - return sphere_fraction(rlow, rhigh, pr); + case GPARAM_ANG1 : + ang1 = val; + break; - case PMODEL_SCGAUSSIAN : - return gaussian_fraction(rlow, rhigh, pr); + case GPARAM_ANG2 : + ang2 = val; + break; - default : - ERROR("No pmodel in volume_fraction!\n"); - return 1.0; - } -} + case GPARAM_R : + R = val; + break; + case GPARAM_WAVELENGTH : + lambda = val; + break; -/* Return the gradient of "fx" wrt parameter 'k' given the current - * status of the crystal. */ -double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) -{ - double glow, ghigh; - double rlow, rhigh, p; - struct image *image = crystal_get_image(cr); - double R = crystal_get_profile_radius(cr); - double gr; + default : + ERROR("Don't understand parameter %i\n", rv[i]); + break; - get_partial(refl, &rlow, &rhigh, &p); + } + } - if ( k == GPARAM_R ) { + cell = rotate_cell_xy(crystal_get_cell_const(cr), ang1, ang2); + crystal_set_cell(cr, cell); - double Rglow, Rghigh; - double D, psph; + crystal_set_profile_radius(cr, R); + crystal_get_image(cr)->lambda = lambda; +} - D = rlow - rhigh; - psph = volume_fraction(rlow, rhigh, R, pmodel); - Rglow = volume_fraction_rgradient(rlow, R, pmodel); - Rghigh = volume_fraction_rgradient(rhigh, R, pmodel); +static double residual_f(const gsl_vector *v, void *pp) +{ + struct rf_priv *pv = pp; + RefList *list; + struct image im; + Crystal *cr; + UnitCell *cell; + double res; + int i; - gr = 4.0*psph/(3.0*D) + (4.0*R/(3.0*D))*(Rglow - Rghigh); - return gr; + for ( i=0; i<v->size; i++ ) { + if ( gsl_vector_get(v, i) > 100.0 ) return GSL_NAN; + } + cr = crystal_copy(pv->cr); + cell = cell_new_from_cell(crystal_get_cell(cr)); + if ( cell == NULL ) return GSL_NAN; + crystal_set_cell(cr, cell); + im = *crystal_get_image(cr); + crystal_set_image(cr, &im); + apply_parameters(v, pv->initial, pv->rv, cr); + + if ( fabs(crystal_get_profile_radius(cr)) > 5e9 ) { + crystal_free(cr); + if ( pv->verbose ) STATUS("radius > 5e9\n"); + return GSL_NAN; } - /* Calculate the gradient of partiality wrt excitation error. */ - glow = partiality_gradient(rlow, R, pmodel, rlow, rhigh); - ghigh = partiality_gradient(rhigh, R, pmodel, rlow, rhigh); + /* Can happen with grid scans and certain --force-radius values */ + if ( fabs(crystal_get_profile_radius(cr)) < 0.0000001e9 ) { + crystal_free(cr); + if ( pv->verbose ) STATUS("radius very small\n"); + return GSL_NAN; + } - if ( k == GPARAM_DIV ) { + if ( im.lambda <= 0.0 ) { + crystal_free(cr); + if ( pv->verbose ) STATUS("lambda < 0\n"); + return GSL_NAN; + } - double D, psph, ds; - signed int hs, ks, ls; + list = copy_reflist(crystal_get_reflections(cr)); + crystal_set_reflections(cr, list); - D = rlow - rhigh; - psph = volume_fraction(rlow, rhigh, R, pmodel); - get_symmetric_indices(refl, &hs, &ks, &ls); - ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); + update_predictions(cr); + calculate_partialities(cr, PMODEL_XSPHERE); - gr = (ds/2.0)*(glow+ghigh) - 4.0*R*psph*ds/(3.0*D*D); - return gr; + res = residual(cr, pv->full, 0, NULL, NULL, 0); - } + cell_free(crystal_get_cell(cr)); + reflist_free(crystal_get_reflections(cr)); + crystal_free(cr); - gr = r_gradient(crystal_get_cell(cr), k, refl, image) * (glow-ghigh); - return gr; + return res; } -static void apply_cell_shift(UnitCell *cell, int k, double shift) +static double get_initial_param(Crystal *cr, enum gparam p) { - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double as, bs, cs; + switch ( p ) { - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - as = modulus(asx, asy, asz); - bs = modulus(bsx, bsy, bsz); - cs = modulus(csx, csy, csz); + case GPARAM_ANG1 : return 0.0; + case GPARAM_ANG2 : return 0.0; + case GPARAM_R : return crystal_get_profile_radius(cr); + case GPARAM_WAVELENGTH : return crystal_get_image(cr)->lambda; - /* Forbid any step which looks too large */ - switch ( k ) - { - case GPARAM_ASX : - case GPARAM_ASY : - case GPARAM_ASZ : - if ( fabs(shift) > 0.1*as ) return; - break; - - case GPARAM_BSX : - case GPARAM_BSY : - case GPARAM_BSZ : - if ( fabs(shift) > 0.1*bs ) return; - break; - - case GPARAM_CSX : - case GPARAM_CSY : - case GPARAM_CSZ : - if ( fabs(shift) > 0.1*cs ) return; - break; - } - - switch ( k ) - { - case GPARAM_ASX : asx += shift; break; - case GPARAM_ASY : asy += shift; break; - case GPARAM_ASZ : asz += shift; break; - case GPARAM_BSX : bsx += shift; break; - case GPARAM_BSY : bsy += shift; break; - case GPARAM_BSZ : bsz += shift; break; - case GPARAM_CSX : csx += shift; break; - case GPARAM_CSY : csy += shift; break; - case GPARAM_CSZ : csz += shift; break; - } - - cell_set_reciprocal(cell, asx, asy, asz, - bsx, bsy, bsz, - csx, csy, csz); + default: return 0.0; + + } } -/* Apply the given shift to the 'k'th parameter of 'image'. */ -static void apply_shift(Crystal *cr, int k, double shift) +static int check_angle_shifts(gsl_vector *cur, const gsl_vector *initial, + enum gparam *rv, struct rf_priv *residual_f_priv) { - double t; - struct image *image = crystal_get_image(cr); - - switch ( k ) { - - case GPARAM_DIV : - if ( shift > 0.1*image->div ) return; - image->div += shift; - if ( image->div < 0.0 ) image->div = 0.0; - break; - - case GPARAM_R : - t = crystal_get_profile_radius(cr); - if ( shift > 0.1*t ) return; - t += shift; - crystal_set_profile_radius(cr, t); - break; - - case GPARAM_ASX : - case GPARAM_ASY : - case GPARAM_ASZ : - case GPARAM_BSX : - case GPARAM_BSY : - case GPARAM_BSZ : - case GPARAM_CSX : - case GPARAM_CSY : - case GPARAM_CSZ : - apply_cell_shift(crystal_get_cell(cr), k, shift); - break; + int i = 0; + double ang = 0.0; - default : - ERROR("No shift defined for parameter %i\n", k); - abort(); + while ( rv[i] != GPARAM_EOL ) { + if ( (rv[i] == GPARAM_ANG1) || (rv[i] == GPARAM_ANG2) ) { + ang += fabs(get_actual_val(cur, initial, rv, i)); + } + rv++; + } + if ( rad2deg(ang) > 5.0 ) { + ERROR("More than 5 degrees total rotation!\n"); + residual_f_priv->verbose = 1; + double res = residual_f(cur, residual_f_priv); + STATUS("residual after rotation = %e\n", res); + residual_f_priv->verbose = 2; + res = residual_f(initial, residual_f_priv); + STATUS("residual before rotation = %e\n", res); + return 1; } + return 0; } -/* Perform one cycle of post refinement on 'image' against 'full' */ -static double pr_iterate(Crystal *cr, const RefList *full, - PartialityModel pmodel, - int *n_filtered, int verbose) +static RefList *reindex_reflections(RefList *input, SymOpList *amb, + SymOpList *sym, int idx) { - gsl_matrix *M; - gsl_vector *v; - gsl_vector *shifts; - int param; + RefList *n; Reflection *refl; RefListIterator *iter; - RefList *reflections; - double max_shift; - int nref = 0; - int num_params = 0; - enum gparam rv[32]; - double G, B; - - if ( n_filtered != NULL ) *n_filtered = 0; - - rv[num_params++] = GPARAM_ASX; - rv[num_params++] = GPARAM_ASY; - rv[num_params++] = GPARAM_ASZ; - rv[num_params++] = GPARAM_BSX; - rv[num_params++] = GPARAM_BSY; - rv[num_params++] = GPARAM_BSZ; - rv[num_params++] = GPARAM_CSX; - rv[num_params++] = GPARAM_CSY; - rv[num_params++] = GPARAM_CSZ; -// rv[num_params++] = GPARAM_R; -// rv[num_params++] = GPARAM_DIV; - - M = gsl_matrix_calloc(num_params, num_params); - v = gsl_vector_calloc(num_params); - - reflections = crystal_get_reflections(cr); - G = crystal_get_osf(cr); - B = crystal_get_Bfac(cr); - - /* Post-refinement terms */ - for ( refl = first_refl(reflections, &iter); + + n = reflist_new(); + + for ( refl = first_refl(input, &iter); refl != NULL; refl = next_refl(refl, iter) ) { - signed int ha, ka, la; - double I_full, delta_I, esd; - double w; - double I_partial; - int k; - double p, L, s; - double fx; - Reflection *match; - double gradients[num_params]; - - if ( get_flag(refl) ) continue; + signed int h, k, l; + Reflection *rn; - get_indices(refl, &ha, &ka, &la); - match = find_refl(full, ha, ka, la); - if ( match == NULL ) continue; - I_full = get_intensity(match); + get_indices(refl, &h, &k, &l); + get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l); + get_asymm(sym, h, k, l, &h, &k, &l); + rn = add_refl(n, h, k, l); + copy_data(rn, refl); + + get_symmetric_indices(rn, &h, &k, &l); + get_equiv(amb, NULL, idx, h, k, l, &h, &k, &l); + set_symmetric_indices(rn, h, k, l); + } - if ( get_redundancy(match) < 2 ) continue; + return n; +} - p = get_partiality(refl); - L = get_lorentz(refl); - I_partial = get_intensity(refl); - esd = get_esd_intensity(refl); - s = resolution(crystal_get_cell(cr), ha, ka, la); - if ( I_partial < 3.0*esd ) continue; +static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx) +{ + signed int h, k, l; + struct rvec na, nb, nc; + struct rvec as, bs, cs; + + cell_get_reciprocal(cell, &as.u, &as.v, &as.w, + &bs.u, &bs.v, &bs.w, + &cs.u, &cs.v, &cs.w); + + get_equiv(amb, NULL, idx, 1, 0, 0, &h, &k, &l); + na.u = as.u*h + bs.u*k + cs.u*l; + na.v = as.v*h + bs.v*k + cs.v*l; + na.w = as.w*h + bs.w*k + cs.w*l; + + get_equiv(amb, NULL, idx, 0, 1, 0, &h, &k, &l); + nb.u = as.u*h + bs.u*k + cs.u*l; + nb.v = as.v*h + bs.v*k + cs.v*l; + nb.w = as.w*h + bs.w*k + cs.w*l; + + get_equiv(amb, NULL, idx, 0, 0, 1, &h, &k, &l); + nc.u = as.u*h + bs.u*k + cs.u*l; + nc.v = as.v*h + bs.v*k + cs.v*l; + nc.w = as.w*h + bs.w*k + cs.w*l; + + cell_set_reciprocal(cell, na.u, na.v, na.w, + nb.u, nb.v, nb.w, + nc.u, nc.v, nc.w); +} - /* Calculate the weight for this reflection */ - w = (s/1e9)*(s/1e9) / (esd*esd); - /* Calculate all gradients for this reflection */ - for ( k=0; k<num_params; k++ ) { - gradients[k] = gradient(cr, rv[k], refl, pmodel); - gradients[k] *= exp(G)*exp(-B*s*s)*I_full/L; - } +void try_reindex(Crystal *crin, const RefList *full, + SymOpList *sym, SymOpList *amb) +{ + RefList *list; + Crystal *cr; + UnitCell *cell; + double residual_original, residual_flipped; + int idx, n; - for ( k=0; k<num_params; k++ ) { + if ( sym == NULL || amb == NULL ) return; - int g; - double v_c, v_curr; + residual_original = residual(crin, full, 0, NULL, NULL, 0); - for ( g=0; g<num_params; g++ ) { + cr = crystal_copy(crin); - double M_c, M_curr; + n = num_equivs(amb, NULL); - /* Matrix is symmetric */ - if ( g > k ) continue; + for ( idx=0; idx<n; idx++ ) { - M_c = w * gradients[g] * gradients[k]; + cell = cell_new_from_cell(crystal_get_cell(crin)); + if ( cell == NULL ) return; + reindex_cell(cell, amb, idx); + crystal_set_cell(cr, cell); - M_curr = gsl_matrix_get(M, k, g); - gsl_matrix_set(M, k, g, M_curr + M_c); - gsl_matrix_set(M, g, k, M_curr + M_c); + list = reindex_reflections(crystal_get_reflections(crin), + amb, sym, idx); + crystal_set_reflections(cr, list); - } + update_predictions(cr); + calculate_partialities(cr, PMODEL_XSPHERE); - fx = exp(G)*p*exp(-B*s*s)*I_full/L; - delta_I = I_partial - fx; - v_c = w * delta_I * gradients[k]; - v_curr = gsl_vector_get(v, k); - gsl_vector_set(v, k, v_curr + v_c); + residual_flipped = residual(cr, full, 0, NULL, NULL, 1); + if ( residual_flipped < residual_original ) { + crystal_set_cell(crin, cell); + crystal_set_reflections(crin, list); + residual_original = residual_flipped; + } else { + cell_free(crystal_get_cell(cr)); + reflist_free(crystal_get_reflections(cr)); } - - nref++; } - if ( verbose ) { - //STATUS("The original equation:\n"); - //show_matrix_eqn(M, v); - STATUS("%i reflections went into the equations.\n", nref); - } - - if ( nref < num_params ) { - crystal_set_user_flag(cr, PRFLAG_FEWREFL); - gsl_matrix_free(M); - gsl_vector_free(v); - return 0.0; - } - - max_shift = 0.0; - shifts = solve_svd(v, M, n_filtered, 0); - if ( shifts != NULL ) { - - for ( param=0; param<num_params; param++ ) { - double shift = gsl_vector_get(shifts, param); - if ( verbose ) STATUS("Shift %i: %e\n", param, shift); - if ( isnan(shift) ) { - //ERROR("NaN shift parameter %i (image ser %i), " - // "%i reflections.\n", rv[param], - // crystal_get_image(cr)->serial, - // nref); - } else { - apply_shift(cr, rv[param], shift); - } - if ( fabs(shift) > max_shift ) max_shift = fabs(shift); - } + crystal_free(cr); +} + + +void write_test_logs(Crystal *crystal, const RefList *full, + signed int cycle, int serial) +{ + FILE *fh; + struct image *image = crystal_get_image(crystal); + char tmp[256]; + char ins[5]; + + snprintf(tmp, 256, "pr-logs/parameters-crystal%i.dat", serial); + + if ( cycle == 0 ) { + fh = fopen(tmp, "w"); } else { - crystal_set_user_flag(cr, PRFLAG_SOLVEFAIL); + fh = fopen(tmp, "a"); + } + + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", tmp); + return; } - gsl_matrix_free(M); - gsl_vector_free(v); - gsl_vector_free(shifts); + if ( cycle == 0 ) { + fprintf(fh, "Image: %s %s\n", + image->filename, get_event_string(image->event)); + } - return max_shift; + if ( cycle >= 0 ) { + snprintf(ins, 4, "%i", cycle); + } else { + ins[0] = 'F'; + ins[1] = '\0'; + } + + fprintf(fh, "%s rlp_size = %e\n", ins, crystal_get_profile_radius(crystal)/1e10); + fprintf(fh, "%s mosaicity = %e\n", ins, crystal_get_mosaicity(crystal)); + fprintf(fh, "%s wavelength = %f A\n", ins, image->lambda*1e10); + fprintf(fh, "%s bandwidth = %f\n", ins, image->bw); + fprintf(fh, "%s my scale factor = %e\n", ins, crystal_get_osf(crystal)); + + double asx, asy, asz, bsx, bsy, bsz, csx, csy, csz; + cell_get_reciprocal(crystal_get_cell(crystal), &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); + fprintf(fh, "%s %f, %f, %f, %f, %f, %f, %f, %f, %f\n", ins, + asx/1e10, bsx/1e10, csx/1e10, + asy/1e10, bsy/1e10, csy/1e10, + asz/1e10, bsz/1e10, csz/1e10); + fclose(fh); } -double residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename) +void write_specgraph(Crystal *crystal, const RefList *full, + signed int cycle, int serial) { - double dev = 0.0; - double G, B; + FILE *fh; + char tmp[256]; Reflection *refl; RefListIterator *iter; - FILE *fh = NULL; - int n_used = 0; + double G = crystal_get_osf(crystal); + double B = crystal_get_Bfac(crystal); + UnitCell *cell; + struct image *image = crystal_get_image(crystal); + char ins[5]; - if ( filename != NULL ) { - fh = fopen(filename, "a"); - if ( fh == NULL ) { - ERROR("Failed to open '%s'\n", filename); - } + snprintf(tmp, 256, "pr-logs/specgraph-crystal%i.dat", serial); + + if ( cycle == 0 ) { + fh = fopen(tmp, "w"); + } else { + fh = fopen(tmp, "a"); } - G = crystal_get_osf(cr); - B = crystal_get_Bfac(cr); + if ( fh == NULL ) { + ERROR("Failed to open '%s'\n", tmp); + return; + } - for ( refl = first_refl(crystal_get_reflections(cr), &iter); + if ( cycle == 0 ) { + fprintf(fh, "Image: %s %s\n", + image->filename, get_event_string(image->event)); + fprintf(fh, "khalf/m 1/d(m) pcalc pobs iteration h k l\n"); + } + + cell = crystal_get_cell(crystal); + + if ( cycle >= 0 ) { + snprintf(ins, 4, "%i", cycle); + } else { + ins[0] = 'F'; + ins[1] = '\0'; + } + + for ( refl = first_refl(crystal_get_reflections(crystal), &iter); refl != NULL; refl = next_refl(refl, iter) ) { - double p, L, s, w; + double corr, Ipart, Ifull, pobs, pcalc; + double res; signed int h, k, l; Reflection *match; - double esd, I_full, I_partial; - double fx, dc; - - if ( free && !get_flag(refl) ) continue; get_indices(refl, &h, &k, &l); + res = resolution(cell, h, k, l); + match = find_refl(full, h, k, l); if ( match == NULL ) continue; - I_full = get_intensity(match); - if ( get_redundancy(match) < 2 ) continue; + /* Don't calculate pobs if reference reflection is weak */ + if ( fabs(get_intensity(match)) / get_esd_intensity(match) < 3.0 ) continue; - p = get_partiality(refl); - L = get_lorentz(refl); - I_partial = get_intensity(refl); - esd = get_esd_intensity(refl); - s = resolution(crystal_get_cell(cr), h, k, l); - - if ( I_partial < 3.0*esd ) continue; + corr = G * exp(B*res*res) * get_lorentz(refl); + Ipart = get_intensity(refl) * corr; + Ifull = get_intensity(match); + pobs = Ipart / Ifull; + pcalc = get_partiality(refl); - fx = exp(G)*p*exp(-B*s*s)*I_full/L; - dc = I_partial - fx; - w = (s/1e9)*(s/1e9)/(esd*esd); - dev += w*dc*dc; - n_used++; + fprintf(fh, "%e %e %f %f %s %4i %4i %4i\n", + get_khalf(refl), 2.0*res, pcalc, pobs, ins, h, k, l); - if ( fh != NULL ) { - fprintf(fh, "%4i %4i %4i %e %e\n", - h, k, l, s, dev); - } } - if ( fh != NULL ) fclose(fh); - - if ( pn_used != NULL ) *pn_used = n_used; - return dev; + fclose(fh); } -static void write_residual_graph(Crystal *cr, const RefList *full) +static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full, + int verbose, int serial, + struct rf_priv *priv) { + gsl_multimin_fminimizer *min; + int n_params = 0; int i; - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - double a; - double step; - double orig_asx; - FILE *fh; - UnitCell *cell; - cell = crystal_get_cell(cr); + /* The parameters to be refined */ + priv->rv[n_params++] = GPARAM_ANG1; + priv->rv[n_params++] = GPARAM_ANG2; + priv->rv[n_params++] = GPARAM_R; + priv->rv[n_params++] = GPARAM_WAVELENGTH; + priv->rv[n_params] = GPARAM_EOL; /* End of list */ + + priv->cr = cr; + priv->full = full; + priv->verbose = verbose; + priv->serial = serial; + + priv->f.f = residual_f; + priv->f.n = n_params; + priv->f.params = priv; + + priv->initial = gsl_vector_alloc(n_params); + priv->vals = gsl_vector_alloc(n_params); + priv->step = gsl_vector_alloc(n_params); + + for ( i=0; i<n_params; i++ ) { + gsl_vector_set(priv->initial, i, get_initial_param(cr, priv->rv[i])); + gsl_vector_set(priv->vals, i, 0.0); + gsl_vector_set(priv->step, i, 1.0); + } - fh = fopen("residual-plot.dat", "w"); + min = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex2, + n_params); + gsl_multimin_fminimizer_set(min, &priv->f, priv->vals, priv->step); - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - a = modulus(asx, asy, asz); - orig_asx = asx; + return min; +} - step = a/100.0; - for ( i=-50; i<=50; i++ ) { +static void write_grid(Crystal *cr, const RefList *full, + signed int cycle, int serial, + const enum gparam par1, const enum gparam par2, + const char *name) +{ + FILE *fh; + char fn[64]; + char ins[5]; + gsl_multimin_fminimizer *min; + struct rf_priv priv; + int idx1, idx2; + int i; - double res; - int n; - asx = orig_asx + (i*step); - cell_set_reciprocal(cell, asx, asy, asz, - bsx, bsy, bsz, - csx, csy, csz); - update_predictions(cr); - calculate_partialities(cr, PMODEL_SCSPHERE); - res = residual(cr, full, 0, &n, NULL); - fprintf(fh, "%i %e %e %i\n", i, asx, res, n); + min = setup_minimiser(cr, full, 0, serial, &priv); + + idx1 = 99; + idx2 = 99; + for ( i=0; i<priv.f.n; i++ ) { + if ( priv.rv[i] == par1 ) idx1 = i; + if ( priv.rv[i] == par2 ) idx2 = i; } + assert(idx1 != 99); + assert(idx2 != 99); - cell_set_reciprocal(cell, orig_asx, asy, asz, - bsx, bsy, bsz, - csx, csy, csz); - update_predictions(cr); - calculate_partialities(cr, PMODEL_SCSPHERE); + if ( cycle >= 0 ) { + snprintf(ins, 4, "%i", cycle); + } else { + ins[0] = 'F'; + ins[1] = '\0'; + } + + snprintf(fn, 63, "pr-logs/grid-crystal%i-cycle%s-%s.dat", + serial, ins, name); + fh = fopen(fn, "w"); + if ( fh != NULL ) { + double v1, v2; + fprintf(fh, "%e %e %e %s\n", + -5.0*get_scale(par1)+get_initial_param(cr, par1), + 5.0*get_scale(par1)+get_initial_param(cr, par1), + get_initial_param(cr, par1), + get_label(par1)); + fprintf(fh, "%e %e %e %s\n", + -5.0*get_scale(par2)+get_initial_param(cr, par2), + 5.0*get_scale(par2)+get_initial_param(cr, par2), + get_initial_param(cr, par2), + get_label(par2)); + for ( v2=-5.0; v2<=5.0; v2+=0.25 ) { + int first=1; + for ( v1=-5.0; v1<=5.0; v1+=0.25 ) { + double res; + gsl_vector_set(min->x, idx1, v1); + gsl_vector_set(min->x, idx2, v2); + res = residual_f(min->x, &priv); + if ( !first ) fprintf(fh, " "); + first = 0; + fprintf(fh, "%e", res); + } + fprintf(fh, "\n"); + } + } fclose(fh); + + gsl_multimin_fminimizer_free(min); + gsl_vector_free(priv.initial); + gsl_vector_free(priv.vals); + gsl_vector_free(priv.step); +} + + +void write_gridscan(Crystal *cr, const RefList *full, + signed int cycle, int serial) +{ + write_grid(cr, full, cycle, serial, + GPARAM_ANG1, GPARAM_ANG2, "ang1-ang2"); + write_grid(cr, full, cycle, serial, + GPARAM_ANG1, GPARAM_WAVELENGTH, "ang1-wave"); + write_grid(cr, full, cycle, serial, + GPARAM_R, GPARAM_WAVELENGTH, "R-wave"); } static void do_pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int verbose) + PartialityModel pmodel, int verbose, int serial, + int cycle, int write_logs, + SymOpList *sym, SymOpList *amb) { - int i, done; - double old_dev; - UnitCell *cell = crystal_get_cell(cr); + gsl_multimin_fminimizer *min; + struct rf_priv priv; + int n_iter = 0; + int status; + int r; + double G; + double residual_start, residual_free_start; + FILE *fh = NULL; + + try_reindex(cr, full, sym, amb); - old_dev = residual(cr, full, 0, 0, NULL); + residual_start = residual(cr, full, 0, NULL, NULL, 0); + residual_free_start = residual(cr, full, 1, NULL, NULL, 0); if ( verbose ) { - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - STATUS("Initial asx = %e\n", asx); - STATUS("PR initial dev = %10.5e, free dev = %10.5e\n", - old_dev, residual(cr, full, 1, NULL, NULL)); - } - - i = 0; - done = 0; - do { + STATUS("\nPR initial: dev = %10.5e, free dev = %10.5e\n", + residual_start, residual_free_start); + } - double dev; + min = setup_minimiser(cr, full, verbose, serial, &priv); - pr_iterate(cr, full, pmodel, NULL, verbose); + if ( verbose ) { + double res = residual_f(min->x, &priv); + double size = gsl_multimin_fminimizer_size(min); + STATUS("At start: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n", + gsl_vector_get(min->x, 0), + gsl_vector_get(min->x, 1), + gsl_vector_get(min->x, 2), + gsl_vector_get(min->x, 3), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), + get_actual_val(min->x, priv.initial, priv.rv, 2), + get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10, + res, size); + } - update_predictions(cr); - calculate_partialities(cr, pmodel); + if ( write_logs ) { - dev = residual(cr, full, 0, 0, NULL); - if ( fabs(dev - old_dev) < dev*0.0001 ) done = 1; + char fn[64]; + + snprintf(fn, 63, "pr-logs/crystal%i-cycle%i.log", serial, cycle); + fh = fopen(fn, "w"); + if ( fh != NULL ) { + fprintf(fh, "iteration RtoReference CCtoReference nref " + "ang1 ang2 radius wavelength"); + double res = residual_f(min->x, &priv); + fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n", + n_iter, res, 0.0, 0, + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), + get_actual_val(min->x, priv.initial, priv.rv, 2), + get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10); + } + + } + + do { + double res; + + n_iter++; + + status = gsl_multimin_fminimizer_iterate(min); + if ( status ) break; + + res = residual_f(min->x, &priv); + if ( isnan(res) ) { + status = GSL_ENOPROG; + break; + } if ( verbose ) { - STATUS("PR iter %2i: dev = %10.5e, free dev = %10.5e\n", - i+1, dev, residual(cr, full, 1, NULL, NULL)); + double res = residual_f(min->x, &priv); + double size = gsl_multimin_fminimizer_size(min); + STATUS("%f %f %f %f ----> %f %f %e %f residual = %e size %f\n", + gsl_vector_get(min->x, 0), + gsl_vector_get(min->x, 1), + gsl_vector_get(min->x, 2), + gsl_vector_get(min->x, 3), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), + get_actual_val(min->x, priv.initial, priv.rv, 2), + get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10, + res, size); } - i++; - old_dev = dev; + if ( fh != NULL ) { + fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n", + n_iter, res, 0.0, 0, + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), + get_actual_val(min->x, priv.initial, priv.rv, 2), + get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10); + } - } while ( i < 30 && !done ); + status = gsl_multimin_test_size(min->size, 0.005); + + } while ( status == GSL_CONTINUE && n_iter < 1000 ); if ( verbose ) { - double asx, asy, asz; - double bsx, bsy, bsz; - double csx, csy, csz; - cell_get_reciprocal(cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); - STATUS("Final asx = %e\n", asx); + STATUS("Done with refinement after %i iter\n", n_iter); + STATUS("status = %i (%s)\n", status, gsl_strerror(status)); } -} + if ( status == GSL_SUCCESS ) { -static struct prdata pr_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel) -{ - int verbose = 0; - struct prdata prdata; + if ( check_angle_shifts(min->x, priv.initial, priv.rv, &priv) ) return; - prdata.refined = 0; - prdata.n_filtered = 0; + if ( verbose ) { - if ( verbose ) { - write_residual_graph(cr, full); + double res = residual_f(min->x, &priv); + double size = gsl_multimin_fminimizer_size(min); + STATUS("At end: %f %f %f %f ----> %f %f %e %f residual = %e size %f\n", + gsl_vector_get(min->x, 0), + gsl_vector_get(min->x, 1), + gsl_vector_get(min->x, 2), + gsl_vector_get(min->x, 3), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), + get_actual_val(min->x, priv.initial, priv.rv, 2), + get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10, + res, size); + + } + + if ( fh != NULL ) { + double res = residual_f(min->x, &priv); + fprintf(fh, "%5i %10.8f %10.8f %5i %10.8f %10.8f %e %e\n", + n_iter, res, 0.0, 0, + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 0)), + rad2deg(get_actual_val(min->x, priv.initial, priv.rv, 1)), + get_actual_val(min->x, priv.initial, priv.rv, 2), + get_actual_val(min->x, priv.initial, priv.rv, 3)*1e10); + fclose(fh); + } + + /* Apply the final shifts */ + apply_parameters(min->x, priv.initial, priv.rv, cr); + update_predictions(cr); + calculate_partialities(cr, PMODEL_XSPHERE); + r = linear_scale(full, crystal_get_reflections(cr), &G, 0); + if ( r == 0 ) { + crystal_set_osf(cr, G); + } + + if ( verbose ) { + + STATUS("After applying final shifts:\n"); + STATUS("PR final: dev = %10.5e, free dev = %10.5e\n", + residual(cr, full, 0, NULL, NULL, 0), + residual(cr, full, 1, NULL, NULL, 0)); + STATUS("Final R = %e m^-1\n", crystal_get_profile_radius(cr)); + + } + + } else { + ERROR("Bad refinement: crystal %i\n", serial); } - do_pr_refine(cr, full, pmodel, verbose); + if ( write_logs ) { + write_gridscan(cr, full, cycle, serial); + write_specgraph(cr, full, cycle, serial); + write_test_logs(cr, full, cycle, serial); + } - if ( crystal_get_user_flag(cr) == 0 ) { - prdata.refined = 1; + if ( crystal_get_profile_radius(cr) > 5e9 ) { + ERROR("Very large radius: crystal %i\n", serial); } - return prdata; + gsl_multimin_fminimizer_free(min); + gsl_vector_free(priv.initial); + gsl_vector_free(priv.vals); + gsl_vector_free(priv.step); } @@ -706,7 +943,13 @@ struct refine_args RefList *full; Crystal *crystal; PartialityModel pmodel; + int serial; struct prdata prdata; + int verbose; + int cycle; + int no_logs; + SymOpList *sym; + SymOpList *amb; }; @@ -724,8 +967,18 @@ static void refine_image(void *task, int id) { struct refine_args *pargs = task; Crystal *cr = pargs->crystal; + int write_logs = 0; + + write_logs = !pargs->no_logs && (pargs->serial % 20 == 0); + pargs->prdata.refined = 0; - pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel); + do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose, + pargs->serial, pargs->cycle, write_logs, + pargs->sym, pargs->amb); + + if ( crystal_get_user_flag(cr) == 0 ) { + pargs->prdata.refined = 1; + } } @@ -738,6 +991,7 @@ static void *get_image(void *vqargs) memcpy(task, &qargs->task_defaults, sizeof(struct refine_args)); task->crystal = qargs->crystals[qargs->n_started]; + task->serial = qargs->n_started; qargs->n_started++; @@ -757,7 +1011,9 @@ static void done_image(void *vqargs, void *task) void refine_all(Crystal **crystals, int n_crystals, - RefList *full, int nthreads, PartialityModel pmodel) + RefList *full, int nthreads, PartialityModel pmodel, + int verbose, int cycle, int no_logs, + SymOpList *sym, SymOpList *amb) { struct refine_args task_defaults; struct queue_args qargs; @@ -766,7 +1022,11 @@ void refine_all(Crystal **crystals, int n_crystals, task_defaults.crystal = NULL; task_defaults.pmodel = pmodel; task_defaults.prdata.refined = 0; - task_defaults.prdata.n_filtered = 0; + task_defaults.verbose = verbose; + task_defaults.cycle = cycle; + task_defaults.no_logs = no_logs; + task_defaults.sym = sym; + task_defaults.amb = amb; qargs.task_defaults = task_defaults; qargs.n_started = 0; diff --git a/src/post-refinement.h b/src/post-refinement.h index 7f395013..e2416a23 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2018 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -41,6 +41,7 @@ #include "utils.h" #include "crystal.h" #include "geometry.h" +#include "symmetry.h" enum prflag @@ -51,19 +52,31 @@ enum prflag PRFLAG_EARLY = 18, PRFLAG_CC = 19, PRFLAG_BIGB = 20, + PRFLAG_SCALEBAD = 21, }; extern const char *str_prflag(enum prflag flag); extern void refine_all(Crystal **crystals, int n_crystals, - RefList *full, int nthreads, PartialityModel pmodel); + RefList *full, int nthreads, PartialityModel pmodel, + int verbose, int cycle, int no_logs, + SymOpList *sym, SymOpList *amb); + +extern void write_gridscan(Crystal *cr, const RefList *full, + int cycle, int serial); + +extern void write_specgraph(Crystal *crystal, const RefList *full, + signed int cycle, int serial); /* Exported so it can be poked by tests/pr_p_gradient_check */ extern double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel); extern double residual(Crystal *cr, const RefList *full, int free, - int *pn_used, const char *filename); + int *pn_used, const char *filename, int complain); + +extern void write_test_logs(Crystal *crystal, const RefList *full, + signed int cycle, int serial); #endif /* POST_REFINEMENT_H */ diff --git a/src/process_image.c b/src/process_image.c index 87089289..dbd4f427 100644 --- a/src/process_image.c +++ b/src/process_image.c @@ -246,7 +246,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, if ( iargs->fix_bandwidth >= 0.0 ) { image.bw = iargs->fix_bandwidth; } else { - image.bw = 0.00000001; + image.bw = 0.0013; } if ( image_feature_count(image.features) < iargs->min_peaks ) { @@ -301,7 +301,7 @@ void process_image(const struct index_args *iargs, struct pattern_args *pargs, /* Integrate! */ time_accounts_set(taccs, TACC_INTEGRATION); sb_shared->pings[cookie]++; - integrate_all_5(&image, iargs->int_meth, PMODEL_SCSPHERE, + integrate_all_5(&image, iargs->int_meth, PMODEL_XSPHERE, iargs->push_res, iargs->ir_inn, iargs->ir_mid, iargs->ir_out, iargs->int_diag, iargs->int_diag_h, diff --git a/src/rejection.c b/src/rejection.c index 10f13f9c..9d41f9b4 100644 --- a/src/rejection.c +++ b/src/rejection.c @@ -137,7 +137,7 @@ static void check_cc(Crystal *cr, RefList *full) pcalc = get_partiality(refl); /* Observed partiality */ - corr = exp(-G) * exp(B*res*res) * get_lorentz(refl); + corr = G * exp(B*res*res) * get_lorentz(refl); Ipart = get_intensity(refl) * corr; pobs = Ipart / get_intensity(match); diff --git a/src/scaling.c b/src/scaling.c index 77a64384..f4c93c66 100644 --- a/src/scaling.c +++ b/src/scaling.c @@ -3,11 +3,11 @@ * * Scaling * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -38,16 +38,14 @@ #include <gsl/gsl_linalg.h> #include <gsl/gsl_eigen.h> #include <gsl/gsl_blas.h> +#include <gsl/gsl_fit.h> #include "merge.h" #include "post-refinement.h" #include "symmetry.h" #include "cell.h" #include "cell-utils.h" - - -/* Maximum number of iterations of NLSq to do for each image per macrocycle. */ -#define MAX_CYCLES (10) +#include "scaling.h" /* Apply the given shift to the 'k'th parameter of 'image'. */ @@ -64,9 +62,9 @@ static void apply_shift(Crystal *cr, int k, double shift) break; case GPARAM_OSF : - t = crystal_get_osf(cr); + t = -log(crystal_get_osf(cr)); t += shift; - crystal_set_osf(cr, t); + crystal_set_osf(cr, exp(-t)); break; default : @@ -78,8 +76,7 @@ static void apply_shift(Crystal *cr, int k, double shift) /* Perform one cycle of scaling of 'cr' against 'full' */ -static double scale_iterate(Crystal *cr, const RefList *full, - PartialityModel pmodel, int *nr) +static double scale_iterate(Crystal *cr, const RefList *full, int *nr) { gsl_matrix *M; gsl_vector *v; @@ -182,7 +179,7 @@ static double scale_iterate(Crystal *cr, const RefList *full, } - fx = G + log(p) - log(L) - B*s*s + log(I_full); + fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full); delta_I = log(I_partial) - fx; v_c = w * delta_I * gradients[k]; v_curr = gsl_vector_get(v, k); @@ -271,7 +268,7 @@ double log_residual(Crystal *cr, const RefList *full, int free, if ( I_full <= 0 ) continue; /* Because log */ if ( p <= 0.0 ) continue; /* Because of log */ - fx = G + log(p) - log(L) - B*s*s + log(I_full); + fx = -log(G) + log(p) - log(L) - B*s*s + log(I_full); dc = log(I_partial) - fx; w = 1.0; dev += w*dc*dc; @@ -290,8 +287,7 @@ double log_residual(Crystal *cr, const RefList *full, int free, } -static void do_scale_refine(Crystal *cr, const RefList *full, - PartialityModel pmodel, int *nr) +static void do_scale_refine(Crystal *cr, const RefList *full, int *nr) { int i, done; double old_dev; @@ -304,7 +300,7 @@ static void do_scale_refine(Crystal *cr, const RefList *full, double dev; - scale_iterate(cr, full, pmodel, nr); + scale_iterate(cr, full, nr); dev = log_residual(cr, full, 0, 0, NULL); if ( fabs(dev - old_dev) < dev*0.01 ) done = 1; @@ -312,7 +308,7 @@ static void do_scale_refine(Crystal *cr, const RefList *full, i++; old_dev = dev; - } while ( i < MAX_CYCLES && !done ); + } while ( i < 10 && !done ); } @@ -320,7 +316,6 @@ struct scale_args { RefList *full; Crystal *crystal; - PartialityModel pmodel; int n_reflections; }; @@ -339,8 +334,7 @@ struct queue_args static void scale_crystal(void *task, int id) { struct scale_args *pargs = task; - do_scale_refine(pargs->crystal, pargs->full, pargs->pmodel, - &pargs->n_reflections); + do_scale_refine(pargs->crystal, pargs->full, &pargs->n_reflections); } @@ -381,14 +375,12 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full, int n = 0; for ( i=0; i<n_crystals; i++ ) { - double r; if ( crystal_get_user_flag(crystals[i]) ) continue; r = log_residual(crystals[i], full, 0, NULL, NULL); if ( isnan(r) ) continue; total += r; n++; - } if ( ninc != NULL ) *ninc = n; @@ -397,8 +389,7 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full, /* Perform iterative scaling, all the way to convergence */ -void scale_all(Crystal **crystals, int n_crystals, int nthreads, - PartialityModel pmodel) +void scale_all(Crystal **crystals, int n_crystals, int nthreads) { struct scale_args task_defaults; struct queue_args qargs; @@ -406,7 +397,6 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, int niter = 0; task_defaults.crystal = NULL; - task_defaults.pmodel = pmodel; qargs.task_defaults = task_defaults; qargs.n_crystals = n_crystals; @@ -422,7 +412,7 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, double bef_res; full = merge_intensities(crystals, n_crystals, nthreads, - pmodel, 2, INFINITY, 0); + 2, INFINITY, 0); old_res = new_res; bef_res = total_log_r(crystals, n_crystals, full, NULL); @@ -456,3 +446,178 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads, ERROR("Too many iterations - giving up!\n"); } } + + +static void scale_crystal_linear(void *task, int id) +{ + struct scale_args *pargs = task; + int r; + double G; + + /* Simple iterative algorithm */ + r = linear_scale(pargs->full, crystal_get_reflections(pargs->crystal), &G, 0); + if ( r == 0 ) { + crystal_set_osf(pargs->crystal, G); + } else { + crystal_set_user_flag(pargs->crystal, PRFLAG_SCALEBAD); + } +} + + +/* Calculates G, by which list2 should be multiplied to fit list1 */ +int linear_scale(const RefList *list1, const RefList *list2, double *G, + int complain_loudly) +{ + const Reflection *refl1; + const Reflection *refl2; + RefListIterator *iter; + int max_n = 256; + int n = 0; + double *x; + double *y; + double *w; + int r; + double cov11; + double sumsq; + int n_esd1 = 0; + int n_esd2 = 0; + int n_ih1 = 0; + int n_ih2 = 0; + int n_nan1 = 0; + int n_nan2 = 0; + int n_inf1 = 0; + int n_inf2 = 0; + int n_part = 0; + int n_nom = 0; + + x = malloc(max_n*sizeof(double)); + w = malloc(max_n*sizeof(double)); + y = malloc(max_n*sizeof(double)); + if ( (x==NULL) || (y==NULL) || (w==NULL) ) { + ERROR("Failed to allocate memory for scaling.\n"); + return 1; + } + + int nb = 0; + for ( refl2 = first_refl_const(list2, &iter); + refl2 != NULL; + refl2 = next_refl_const(refl2, iter) ) + { + signed int h, k, l; + double Ih1, Ih2; + double esd1, esd2; + nb++; + + get_indices(refl2, &h, &k, &l); + refl1 = find_refl(list1, h, k, l); + if ( refl1 == NULL ) { + n_nom++; + continue; + } + + Ih1 = get_intensity(refl1); + Ih2 = get_intensity(refl2); + esd1 = get_esd_intensity(refl1); + esd2 = get_esd_intensity(refl2); + + /* Problem cases in approximate descending order of severity */ + if ( isnan(Ih1) ) { n_nan1++; continue; } + if ( isinf(Ih1) ) { n_inf1++; continue; } + if ( isnan(Ih2) ) { n_nan2++; continue; } + if ( isinf(Ih2) ) { n_inf2++; continue; } + if ( get_partiality(refl2) < 0.3 ) { n_part++; continue; } + //if ( Ih1 <= 0.0 ) { n_ih1++; continue; } + if ( Ih2 <= 0.0 ) { n_ih2++; continue; } + //if ( Ih1 <= 3.0*esd1 ) { n_esd1++; continue; } + if ( Ih2 <= 3.0*esd2 ) { n_esd2++; continue; } + //if ( get_redundancy(refl1) < 2 ) continue; + + if ( n == max_n ) { + max_n *= 2; + x = realloc(x, max_n*sizeof(double)); + y = realloc(y, max_n*sizeof(double)); + w = realloc(w, max_n*sizeof(double)); + if ( (x==NULL) || (y==NULL) || (w==NULL) ) { + ERROR("Failed to allocate memory for scaling.\n"); + return 1; + } + } + + x[n] = Ih2 / get_partiality(refl2); + y[n] = Ih1; + w[n] = get_partiality(refl2); + n++; + + } + + if ( n < 2 ) { + if ( complain_loudly ) { + ERROR("Not enough reflections for scaling (had %i, but %i remain)\n", nb, n); + if ( n_esd1 ) ERROR("%i reference reflection esd\n", n_esd1); + if ( n_esd2 ) ERROR("%i subject reflection esd\n", n_esd2); + if ( n_ih1 ) ERROR("%i reference reflection intensity\n", n_ih1); + if ( n_ih2 ) ERROR("%i subject reflection intensity\n", n_ih2); + if ( n_nan1 ) ERROR("%i reference reflection nan\n", n_nan1); + if ( n_nan2 ) ERROR("%i subject reflection nan\n", n_nan2); + if ( n_inf1 ) ERROR("%i reference reflection inf\n", n_inf1); + if ( n_inf2 ) ERROR("%i subject reflection inf\n", n_inf2); + if ( n_part ) ERROR("%i subject reflection partiality\n", n_part); + if ( n_nom ) ERROR("%i no match in reference list\n", n_nom); + } + *G = 1.0; + return 1; + } + + r = gsl_fit_wmul(x, 1, w, 1, y, 1, n, G, &cov11, &sumsq); + + if ( r ) { + ERROR("Scaling failed.\n"); + return 1; + } + + if ( isnan(*G) ) { + + if ( complain_loudly ) { + ERROR("Scaling gave NaN (%i pairs)\n", n); + if ( n < 10 ) { + int i; + for ( i=0; i<n; i++ ) { + STATUS("%i %e %e %e\n", i, x[i], y[i], w[n]); + } + } + } + + *G = 1.0; + return 1; + } + + free(x); + free(y); + free(w); + + return 0; +} + + +void scale_all_to_reference(Crystal **crystals, int n_crystals, + RefList *reference, int nthreads) +{ + struct scale_args task_defaults; + struct queue_args qargs; + + task_defaults.crystal = NULL; + + qargs.task_defaults = task_defaults; + qargs.n_crystals = n_crystals; + qargs.crystals = crystals; + + /* Don't have threads which are doing nothing */ + if ( n_crystals < nthreads ) nthreads = n_crystals; + + qargs.task_defaults.full = reference; + qargs.n_started = 0; + qargs.n_done = 0; + qargs.n_reflections = 0; + run_threads(nthreads, scale_crystal_linear, get_crystal, done_crystal, + &qargs, n_crystals, 0, 0, 0); +} diff --git a/src/scaling.h b/src/scaling.h index 534043d9..514bee03 100644 --- a/src/scaling.h +++ b/src/scaling.h @@ -3,11 +3,11 @@ * * Scaling * - * Copyright © 2012-2015 Deutsches Elektronen-Synchrotron DESY, + * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY, * a research centre of the Helmholtz Association. * * Authors: - * 2010-2015 Thomas White <taw@physics.org> + * 2010-2017 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -41,7 +41,12 @@ extern double log_residual(Crystal *cr, const RefList *full, int free, int *pn_used, const char *filename); -extern void scale_all(Crystal **crystals, int n_crystals, int nthreads, - PartialityModel pmodel); +extern int linear_scale(const RefList *list1, const RefList *list2, double *G, + int complain_loudly); + +extern void scale_all(Crystal **crystals, int n_crystals, int nthreads); + +extern void scale_all_to_reference(Crystal **crystals, int n_crystals, + RefList *reference, int nthreads); #endif /* SCALING_H */ |