diff options
author | Thomas White <taw@physics.org> | 2013-05-02 16:43:22 +0200 |
---|---|---|
committer | Thomas White <taw@physics.org> | 2013-05-02 16:43:22 +0200 |
commit | a4bdd1335ffc98e885126af5468b0696396479cf (patch) | |
tree | 43fe08d1ccb8f22acfd962f6302713e9ff2ddf6f /src | |
parent | c5bcfa0ae036145f23c344796f8b4efcee69429a (diff) | |
parent | 6b18fa18fbdc1f272e65f81d04b22d973acafa1e (diff) |
Merge branch 'tom/partials'
Diffstat (limited to 'src')
-rw-r--r-- | src/hrs-scaling.c | 138 | ||||
-rw-r--r-- | src/hrs-scaling.h | 10 | ||||
-rw-r--r-- | src/partial_sim.c | 88 | ||||
-rw-r--r-- | src/partialator.c | 59 | ||||
-rw-r--r-- | src/post-refinement.c | 114 | ||||
-rw-r--r-- | src/post-refinement.h | 6 |
6 files changed, 293 insertions, 122 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 0498a532..c249ef46 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -3,11 +3,11 @@ * * Intensity scaling using generalised HRS target function * - * 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 <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -60,6 +60,7 @@ struct scale_queue_args Crystal **crystals; int n_started; double max_shift; + PartialityModel pmodel; }; @@ -68,6 +69,7 @@ struct scale_worker_args Crystal *crystal; double shift; RefList *reference; + PartialityModel pmodel; }; @@ -78,6 +80,7 @@ static void *create_scale_job(void *vqargs) wargs = malloc(sizeof(struct scale_worker_args)); wargs->reference = qargs->reference; + wargs->pmodel = qargs->pmodel; wargs->crystal = qargs->crystals[qargs->n_started++]; @@ -94,8 +97,8 @@ static void run_scale_job(void *vwargs, int cookie) RefListIterator *iter; double num = 0.0; double den = 0.0; - double corr; - const double osf = crystal_get_osf(cr); + double g; + const double G = crystal_get_osf(cr); if ( crystal_get_user_flag(cr) ) { wargs->shift = 0.0; @@ -109,6 +112,7 @@ static void run_scale_job(void *vwargs, int cookie) signed int h, k, l; double Ih, Ihl, esd; Reflection *r; + double corr; if ( !get_scalable(refl) ) continue; @@ -125,24 +129,41 @@ static void run_scale_job(void *vwargs, int cookie) Ih = get_intensity(r); } - /* FIXME: This isn't correct */ - Ihl = get_intensity(refl) / get_partiality(refl); - esd = get_esd_intensity(refl) / get_partiality(refl); + /* If you change this, be sure to also change + * run_merge_job() and run_esd_job(). */ + switch ( wargs->pmodel ) { - num += Ih * (Ihl/osf) / pow(esd/osf, 2.0); - den += pow(Ih, 2.0)/pow(esd/osf, 2.0); + case PMODEL_UNITY : + corr = 1.0; + break; + + case PMODEL_SPHERE : + corr = get_partiality(refl) * get_lorentz(refl); + break; + + default : + ERROR("Unrecognised partiality model!\n"); + abort(); + break; + + } + + Ihl = get_intensity(refl) / corr; + esd = get_esd_intensity(refl) / corr; + + num += Ih * (Ihl/G) / pow(esd/G, 2.0); + den += pow(Ih, 2.0)/pow(esd/G, 2.0); } //num += image->osf / pow(SCALING_RESTRAINT, 2.0); //den += pow(image->osf, 2.0)/pow(SCALING_RESTRAINT, 2.0); - corr = num / den; - if ( !isnan(corr) && !isinf(corr) ) { - crystal_set_osf(cr, osf*corr); + g = num / den; + if ( !isnan(g) && !isinf(g) ) { + crystal_set_osf(cr, g*G); } - wargs->shift = fabs(corr-1.0); - + wargs->shift = fabs(g-1.0); } @@ -157,7 +178,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs) static double iterate_scale(Crystal **crystals, int n, RefList *reference, - int n_threads) + int n_threads, PartialityModel pmodel) { struct scale_queue_args qargs; @@ -167,6 +188,7 @@ static double iterate_scale(Crystal **crystals, int n, RefList *reference, qargs.n_started = 0; qargs.crystals = crystals; qargs.max_shift = 0.0; + qargs.pmodel = pmodel; run_threads(n_threads, run_scale_job, create_scale_job, finalise_scale_job, &qargs, n, 0, 0, 0); @@ -181,6 +203,7 @@ struct merge_queue_args pthread_mutex_t full_lock; Crystal **crystals; int n_started; + PartialityModel pmodel; }; @@ -189,6 +212,7 @@ struct merge_worker_args Crystal *crystal; RefList *full; pthread_mutex_t *full_lock; + PartialityModel pmodel; }; @@ -200,6 +224,7 @@ static void *create_merge_job(void *vqargs) wargs = malloc(sizeof(struct merge_worker_args)); wargs->full = qargs->full; wargs->full_lock = &qargs->full_lock; + wargs->pmodel = qargs->pmodel; wargs->crystal = qargs->crystals[qargs->n_started++]; @@ -228,7 +253,7 @@ static void run_merge_job(void *vwargs, int cookie) signed int h, k, l; double num, den; int red; - double Ihl, esd, pcalc; + double Ihl, esd, corr; if ( !get_scalable(refl) ) continue; @@ -251,9 +276,27 @@ static void run_merge_job(void *vwargs, int cookie) red = get_redundancy(f); } - pcalc = get_partiality(refl); - Ihl = get_intensity(refl) / pcalc; - esd = get_esd_intensity(refl) / pcalc; + /* If you change this, be sure to also change + * run_scale_job() and run_esd_job(). */ + switch ( wargs->pmodel ) { + + case PMODEL_UNITY : + corr = 1.0; + break; + + case PMODEL_SPHERE : + corr = get_partiality(refl) * get_lorentz(refl); + break; + + default : + ERROR("Unrecognised partiality model!\n"); + abort(); + break; + + } + + Ihl = get_intensity(refl) / corr; + esd = get_esd_intensity(refl) / corr; num += (Ihl/G) / pow(esd/G, 2.0); den += 1.0 / pow(esd/G, 2.0); @@ -273,7 +316,8 @@ static void finalise_merge_job(void *vqargs, void *vwargs) } -static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads) +static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads, + PartialityModel pmodel) { RefList *full; struct merge_queue_args qargs; @@ -285,6 +329,7 @@ static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads) qargs.full = full; qargs.n_started = 0; qargs.crystals = crystals; + qargs.pmodel = pmodel; pthread_mutex_init(&qargs.full_lock, NULL); run_threads(n_threads, run_merge_job, create_merge_job, @@ -311,6 +356,7 @@ struct esd_queue_args RefList *full; Crystal **crystals; int n_started; + PartialityModel pmodel; }; @@ -318,6 +364,7 @@ struct esd_worker_args { Crystal *crystal; RefList *full; + PartialityModel pmodel; }; @@ -328,6 +375,7 @@ static void *create_esd_job(void *vqargs) wargs = malloc(sizeof(struct esd_worker_args)); wargs->full = qargs->full; + wargs->pmodel = qargs->pmodel; wargs->crystal = qargs->crystals[qargs->n_started++]; @@ -355,7 +403,7 @@ static void run_esd_job(void *vwargs, int cookie) Reflection *f; signed int h, k, l; double num; - double Ihl, Ih; + double Ihl, Ih, corr; if ( !get_scalable(refl) ) continue; @@ -367,8 +415,27 @@ static void run_esd_job(void *vwargs, int cookie) num = get_temp1(f); + /* If you change this, be sure to also change + * run_scale_job() and run_merge_job(). */ + switch ( wargs->pmodel ) { + + case PMODEL_UNITY : + corr = 1.0; + break; + + case PMODEL_SPHERE : + corr = get_partiality(refl) * get_lorentz(refl); + break;; + + default : + ERROR("Unrecognised partiality model!\n"); + abort(); + break; + + } + Ih = get_intensity(f); - Ihl = get_intensity(refl) / (get_partiality(refl) * G); + Ihl = get_intensity(refl) / (corr * G); num += pow(Ihl - Ih, 2.0); @@ -385,7 +452,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs) static void calculate_esds(Crystal **crystals, int n, RefList *full, - int n_threads, int min_red) + int n_threads, int min_red, PartialityModel pmodel) { struct esd_queue_args qargs; Reflection *refl; @@ -394,6 +461,7 @@ static void calculate_esds(Crystal **crystals, int n, RefList *full, qargs.full = full; qargs.n_started = 0; qargs.crystals = crystals; + qargs.pmodel = pmodel; for ( refl = first_refl(full, &iter); refl != NULL; @@ -425,24 +493,25 @@ static void calculate_esds(Crystal **crystals, int n, RefList *full, /* Scale the stack of images */ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, - int n_threads, int noscale) + int n_threads, int noscale, PartialityModel pmodel, + int min_redundancy) { int i; double max_corr; RefList *full = NULL; - const int min_redundancy = 3; for ( i=0; i<n; i++ ) crystal_set_osf(crystals[i], 1.0); if ( noscale ) { - full = lsq_intensities(crystals, n, n_threads); - calculate_esds(crystals, n, full, n_threads, min_redundancy); + full = lsq_intensities(crystals, n, n_threads, pmodel); + calculate_esds(crystals, n, full, n_threads, min_redundancy, + pmodel); return full; } /* No reference -> create an initial list to refine against */ if ( gref == NULL ) { - full = lsq_intensities(crystals, n, n_threads); + full = lsq_intensities(crystals, n, n_threads, pmodel); } /* Iterate */ @@ -458,27 +527,26 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, reference = full; } - max_corr = iterate_scale(crystals, n, reference, n_threads); + max_corr = iterate_scale(crystals, n, reference, n_threads, + pmodel); //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(crystals, n, n_threads); + full = lsq_intensities(crystals, n, n_threads, pmodel); } - //show_scale_factors(images, n); - i++; } while ( (max_corr > 0.01) && (i < MAX_CYCLES) ); if ( gref != NULL ) { - full = lsq_intensities(crystals, n, n_threads); + full = lsq_intensities(crystals, n, n_threads, pmodel); } /* else we already did it */ - calculate_esds(crystals, n, full, n_threads, min_redundancy); + calculate_esds(crystals, n, full, n_threads, min_redundancy, pmodel); return full; } diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h index 80940347..0979bb98 100644 --- a/src/hrs-scaling.h +++ b/src/hrs-scaling.h @@ -3,11 +3,11 @@ * * Intensity scaling using generalised HRS target function * - * 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 <taw@physics.org> + * 2010-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -37,10 +37,12 @@ #include "crystal.h" #include "reflist.h" +#include "geometry.h" extern RefList *scale_intensities(Crystal **crystals, int n, RefList *reference, int n_threads, - int noscale); + int noscale, PartialityModel pmodel, + int min_redundancy); #endif /* HRS_SCALING_H */ diff --git a/src/partial_sim.c b/src/partial_sim.c index 7c93da1b..87d30f39 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -3,11 +3,11 @@ * * Generate partials for testing scaling * - * 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: - * 2011-2012 Thomas White <taw@physics.org> + * 2011-2013 Thomas White <taw@physics.org> * * This file is part of CrystFEL. * @@ -88,7 +88,8 @@ static void calculate_partials(Crystal *cr, int random_intensities, pthread_mutex_t *full_lock, unsigned long int *n_ref, double *p_hist, - double *p_max, double max_q) + double *p_max, double max_q, double full_stddev, + double noise_stddev) { Reflection *refl; RefListIterator *iter; @@ -100,12 +101,13 @@ static void calculate_partials(Crystal *cr, { signed int h, k, l; Reflection *rfull; - double p, Ip, If; + double L, p, Ip, If; int bin; get_indices(refl, &h, &k, &l); get_asymm(sym, h, k, l, &h, &k, &l); p = get_partiality(refl); + L = get_lorentz(refl); pthread_mutex_lock(full_lock); rfull = find_refl(full, h, k, l); @@ -120,7 +122,7 @@ static void calculate_partials(Crystal *cr, * thing under lock. */ pthread_mutex_lock(full_lock); rfull = add_refl(full, h, k, l); - If = fabs(gaussian_noise(0.0, 1000.0)); + If = fabs(gaussian_noise(0.0, full_stddev)); set_intensity(rfull, If); set_redundancy(rfull, 1); pthread_mutex_unlock(full_lock); @@ -139,7 +141,7 @@ static void calculate_partials(Crystal *cr, } } - Ip = crystal_get_osf(cr) * p * If; + Ip = crystal_get_osf(cr) * L * p * If; res = resolution(crystal_get_cell(cr), h, k, l); bin = NBINS*2.0*res/max_q; @@ -152,10 +154,10 @@ static void calculate_partials(Crystal *cr, res, bin, p); } - Ip = gaussian_noise(Ip, 100.0); + Ip = gaussian_noise(Ip, noise_stddev); set_intensity(refl, Ip); - set_esd_intensity(refl, 100.0); + set_esd_intensity(refl, noise_stddev); } } @@ -183,6 +185,10 @@ static void show_help(const char *s) " -c, --cnoise=<val> Add random noise, with a flat distribution, to the\n" " reciprocal lattice vector components given in the\n" " stream, with maximum error +/- <val> percent.\n" +" --osf-stddev=<val> Set the standard deviation of the scaling factors.\n" +" --full-stddev=<val> Set the standard deviation of the randomly\n" +" generated full intensities, if not using -i.\n" +" --noise-stddev=<val> Set the standard deviation of the noise.\n" "\n" ); } @@ -201,6 +207,9 @@ struct queue_args int random_intensities; UnitCell *cell; double cnoise; + double osf_stddev; + double full_stddev; + double noise_stddev; struct image *template_image; double max_q; @@ -254,6 +263,7 @@ static void run_job(void *vwargs, int cookie) int i; Crystal *cr; RefList *reflections; + double osf; cr = crystal_new(); if ( cr == NULL ) { @@ -263,7 +273,10 @@ static void run_job(void *vwargs, int cookie) wargs->crystal = cr; crystal_set_image(cr, &wargs->image); - crystal_set_osf(cr, gaussian_noise(1.0, 0.3)); + do { + osf = gaussian_noise(1.0, qargs->osf_stddev); + } while ( osf <= 0.0 ); + crystal_set_osf(cr, osf); crystal_set_profile_radius(cr, wargs->image.beam->profile_radius); /* Set up a random orientation */ @@ -284,7 +297,8 @@ static void run_job(void *vwargs, int cookie) qargs->sym, qargs->random_intensities, &qargs->full_lock, wargs->n_ref, wargs->p_hist, wargs->p_max, - qargs->max_q); + qargs->max_q, qargs->full_stddev, + qargs->noise_stddev); /* Give a slightly incorrect cell in the stream */ mess_up_cell(cr, qargs->cnoise); @@ -343,6 +357,9 @@ int main(int argc, char *argv[]) int i; FILE *fh; char *phist_file = NULL; + double osf_stddev = 2.0; + double full_stddev = 1000.0; + double noise_stddev = 20.0; /* Long options */ const struct option longopts[] = { @@ -354,8 +371,13 @@ int main(int argc, char *argv[]) {"geometry", 1, NULL, 'g'}, {"symmetry", 1, NULL, 'y'}, {"save-random", 1, NULL, 'r'}, - {"pgraph", 1, NULL, 2}, {"cnoise", 1, NULL, 'c'}, + + {"pgraph", 1, NULL, 2}, + {"osf-stddev", 1, NULL, 3}, + {"full-stddev", 1, NULL, 4}, + {"noise-stddev", 1, NULL, 5}, + {0, 0, NULL, 0} }; @@ -417,6 +439,45 @@ int main(int argc, char *argv[]) phist_file = strdup(optarg); break; + case 3 : + osf_stddev = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid OSF standard deviation.\n"); + return 1; + } + if ( osf_stddev < 0.0 ) { + ERROR("Invalid OSF standard deviation."); + ERROR(" (must be positive).\n"); + return 1; + } + break; + + case 4 : + full_stddev = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid full standard deviation.\n"); + return 1; + } + if ( full_stddev < 0.0 ) { + ERROR("Invalid full standard deviation."); + ERROR(" (must be positive).\n"); + return 1; + } + break; + + case 5 : + noise_stddev = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid noise standard deviation.\n"); + return 1; + } + if ( noise_stddev < 0.0 ) { + ERROR("Invalid noise standard deviation."); + ERROR(" (must be positive).\n"); + return 1; + } + break; + case 0 : break; @@ -548,6 +609,9 @@ int main(int argc, char *argv[]) qargs.template_image = ℑ qargs.stream = stream; qargs.cnoise = cnoise; + qargs.osf_stddev = osf_stddev; + qargs.full_stddev = full_stddev; + qargs.noise_stddev = noise_stddev; qargs.max_q = largest_q(&image); for ( i=0; i<NBINS; i++ ) { diff --git a/src/partialator.c b/src/partialator.c index a4af3c18..cfa3a6ba 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -79,6 +79,9 @@ static void show_help(const char *s) " -r, --reference=<file> Refine images against reflections in <file>,\n" " instead of taking the mean of the intensity\n" " estimates.\n" +" -m, --model=<model> Specify partiality model.\n" +" --min-measurements=<n> Require at least <n> measurements before a\n" +" reflection appears in the output. Default: 2\n" "\n" " -j <n> Run <n> analyses in parallel.\n"); } @@ -88,6 +91,7 @@ struct refine_args { RefList *full; Crystal *crystal; + PartialityModel pmodel; }; @@ -106,7 +110,7 @@ static void refine_image(void *task, int id) struct refine_args *pargs = task; Crystal *cr = pargs->crystal; - pr_refine(cr, pargs->full); + pr_refine(cr, pargs->full, pargs->pmodel); } @@ -139,13 +143,18 @@ static void done_image(void *vqargs, void *task) static void refine_all(Crystal **crystals, int n_crystals, struct detector *det, - RefList *full, int nthreads) + RefList *full, int nthreads, PartialityModel pmodel) { struct refine_args task_defaults; struct queue_args qargs; + /* If the partiality model is "p=1", this refinement is really, really + * easy... */ + if ( pmodel == PMODEL_UNITY ) return; + task_defaults.full = full; task_defaults.crystal = NULL; + task_defaults.pmodel = pmodel; qargs.task_defaults = task_defaults; qargs.n_started = 0; @@ -314,6 +323,10 @@ int main(int argc, char *argv[]) int noscale = 0; Stream *st; Crystal **crystals; + char *pmodel_str = NULL; + PartialityModel pmodel = PMODEL_SPHERE; + int min_measurements = 2; + char *rval; /* Long options */ const struct option longopts[] = { @@ -326,6 +339,9 @@ int main(int argc, char *argv[]) {"iterations", 1, NULL, 'n'}, {"no-scale", 0, &noscale, 1}, {"reference", 1, NULL, 'r'}, + {"model", 1, NULL, 'm'}, + {"min-measurements", 1, NULL, 2}, + {0, 0, NULL, 0} }; @@ -336,7 +352,7 @@ int main(int argc, char *argv[]) } /* Short options */ - while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:r:j:", + while ((c = getopt_long(argc, argv, "hi:o:g:b:y:n:r:j:m:", longopts, NULL)) != -1) { @@ -370,6 +386,10 @@ int main(int argc, char *argv[]) n_iter = atoi(optarg); break; + case 'm' : + pmodel_str = strdup(optarg); + break; + case 'b' : beam = get_beam_parameters(optarg); if ( beam == NULL ) { @@ -383,6 +403,15 @@ int main(int argc, char *argv[]) reference_file = strdup(optarg); break; + case 2 : + errno = 0; + min_measurements = strtod(optarg, &rval); + if ( *rval != '\0' ) { + ERROR("Invalid value for --min-measurements.\n"); + return 1; + } + break; + case 0 : break; @@ -434,6 +463,17 @@ int main(int argc, char *argv[]) return 1; } + if ( pmodel_str != NULL ) { + if ( strcmp(pmodel_str, "sphere") == 0 ) { + pmodel = PMODEL_SPHERE; + } else if ( strcmp(pmodel_str, "unity") == 0 ) { + pmodel = PMODEL_UNITY; + } else { + ERROR("Unknown partiality model '%s'.\n", pmodel_str); + return 1; + } + } + if ( reference_file != NULL ) { RefList *list; @@ -528,6 +568,7 @@ int main(int argc, char *argv[]) display_progress(n_images, n_crystals); } while ( 1 ); + fprintf(stderr, "\n"); close_stream(st); @@ -544,18 +585,19 @@ int main(int argc, char *argv[]) crystal_set_image(cryst, &images[i]); /* Now it's safe to do the following */ - update_partialities(cryst); + update_partialities(cryst, pmodel); as = crystal_get_reflections(cryst); nobs += select_scalable_reflections(as, reference); } } + STATUS("%i scalable observations.\n", nobs); /* Make initial estimates */ - STATUS("\nPerforming initial scaling.\n"); + STATUS("Performing initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); full = scale_intensities(crystals, n_crystals, reference, - nthreads, noscale); + nthreads, noscale, pmodel, min_measurements); sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf", infile, cmdline); @@ -578,7 +620,7 @@ int main(int argc, char *argv[]) /* Refine the geometry of all patterns to get the best fit */ select_reflections_for_refinement(crystals, n_crystals, comp, have_reference); - refine_all(crystals, n_crystals, det, comp, nthreads); + refine_all(crystals, n_crystals, det, comp, nthreads, pmodel); nobs = 0; for ( j=0; j<n_crystals; j++ ) { @@ -593,7 +635,8 @@ int main(int argc, char *argv[]) /* Re-estimate all the full intensities */ reflist_free(full); full = scale_intensities(crystals, n_crystals, - reference, nthreads, noscale); + reference, nthreads, noscale, pmodel, + min_measurements); sr_iteration(sr, i+1, crystals, n_crystals, full); diff --git a/src/post-refinement.c b/src/post-refinement.c index 1439b148..9e4649a2 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -81,27 +81,27 @@ static double partiality_rgradient(double r, double profile_radius) dpdq = 6.0*(q-pow(q, 2.0)); /* dq/drad */ - dqdrad = 0.5 * (1.0 - r * pow(profile_radius, -2.0)); + dqdrad = -0.5 * r * pow(profile_radius, -2.0); return dpdq * dqdrad; } /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(Crystal *cr, int k, Reflection *refl) +double gradient(Crystal *cr, int k, Reflection *refl, PartialityModel pmodel) { - double ds, azix, aziy; - double ttlow, tthigh, tt; - double nom, den; - double g; + double ds, azi; + double glow, ghigh; double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; double xl, yl, zl; signed int hs, ks, ls; - double r1, r2, p; + double rlow, rhigh, p; int clamp_low, clamp_high; - double klow, khigh; + double philow, phihigh, phi; + double khigh, klow; + double tl, cet, cez; double gr; struct image *image = crystal_get_image(cr); double r = crystal_get_profile_radius(cr); @@ -116,33 +116,40 @@ double gradient(Crystal *cr, int k, Reflection *refl) zl = hs*asz + ks*bsz + ls*csz; ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); - get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); + get_partial(refl, &rlow, &rhigh, &p, &clamp_low, &clamp_high); + /* "low" gives the largest Ewald sphere (wavelength short => k large) + * "high" gives the smallest Ewald sphere (wavelength long => k small) + */ klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); khigh = 1.0/(image->lambda + image->lambda*image->bw/2.0); - ttlow = angle_between(0.0, 0.0, 1.0, xl, yl, zl+klow); - tthigh = angle_between(0.0, 0.0, 1.0, xl, yl, zl+khigh); - if ( (clamp_low == 0) && (clamp_high == 0) ) { - tt = (ttlow+tthigh)/2.0; - } else if ( clamp_high == 0 ) { - tt = tthigh + image->div; - } else if ( clamp_low == 0 ) { - tt = ttlow - image->div; - } else { - tt = 0.0; - /* Gradient should come out as zero in this case */ - } - azix = angle_between(1.0, 0.0, 0.0, xl, yl, 0.0); - aziy = angle_between(0.0, 1.0, 0.0, xl, yl, 0.0); + tl = sqrt(xl*xl + yl*yl); + ds = modulus(xl, yl, zl); + + cet = -sin(image->div/2.0) * klow; + cez = -cos(image->div/2.0) * klow; + philow = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0); + + cet = -sin(image->div/2.0) * khigh; + cez = -cos(image->div/2.0) * khigh; + phihigh = M_PI_2 - angle_between_2d(tl-cet, zl-cez, 1.0, 0.0); + + /* Approximation: philow and phihigh are very similar */ + phi = (philow + phihigh) / 2.0; + + azi = atan2(yl, xl); /* Calculate the gradient of partiality wrt excitation error. */ - g = 0.0; if ( clamp_low == 0 ) { - g -= partiality_gradient(r1, r); + glow = partiality_gradient(rlow, r); + } else { + glow = 0.0; } if ( clamp_high == 0 ) { - g += partiality_gradient(r2, r); + ghigh = partiality_gradient(rhigh, r); + } else { + ghigh = 0.0; } /* For many gradients, just multiply the above number by the gradient @@ -150,57 +157,41 @@ double gradient(Crystal *cr, int k, Reflection *refl) switch ( k ) { case REF_DIV : - gr = 0.0; - if ( clamp_low == 0 ) { - nom = sqrt(2.0) * ds * sin(image->div/2.0); - den = sqrt(1.0 - cos(image->div/2.0)); - gr -= (nom/den) * g; - } - if ( clamp_high == 0 ) { - nom = sqrt(2.0) * ds * sin(image->div/2.0); - den = sqrt(1.0 - cos(image->div/2.0)); - gr += (nom/den) * g; - } - if ( isnan(gr) ) gr = 0.0; /* FIXME: This isn't true (?) */ - return gr / 4.0; /* FIXME: Shameless fudge factor */ + /* Small angle approximation */ + return (ds*glow + ds*ghigh) / 2.0; case REF_R : - g = 0.0; - if ( clamp_low == 0 ) { - g += partiality_rgradient(r1, r); - } - if ( clamp_high == 0 ) { - g += partiality_rgradient(r2, r); - } - return g; + gr = partiality_rgradient(rlow, r); + gr -= partiality_rgradient(rhigh, r); + return gr; /* Cell parameters and orientation */ case REF_ASX : - return hs * sin(tt) * cos(azix) * g; + return hs * sin(phi) * cos(azi) * (ghigh-glow); case REF_BSX : - return ks * sin(tt) * cos(azix) * g; + return ks * sin(phi) * cos(azi) * (ghigh-glow); case REF_CSX : - return ls * sin(tt) * cos(azix) * g; + return ls * sin(phi) * cos(azi) * (ghigh-glow); case REF_ASY : - return hs * sin(tt) * cos(aziy) * g; + return hs * sin(phi) * sin(azi) * (ghigh-glow); case REF_BSY : - return ks * sin(tt) * cos(aziy) * g; + return ks * sin(phi) * sin(azi) * (ghigh-glow); case REF_CSY : - return ls * sin(tt) * cos(aziy) * g; + return ls * sin(phi) * sin(azi) * (ghigh-glow); case REF_ASZ : - return hs * cos(tt) * g; + return hs * cos(phi) * (ghigh-glow); case REF_BSZ : - return ks * cos(tt) * g; + return ks * cos(phi) * (ghigh-glow); case REF_CSZ : - return ls * cos(tt) * g; + return ls * cos(phi) * (ghigh-glow); } @@ -364,7 +355,8 @@ 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(Crystal *cr, const RefList *full) +static double pr_iterate(Crystal *cr, const RefList *full, + PartialityModel pmodel) { gsl_matrix *M; gsl_vector *v; @@ -420,7 +412,7 @@ static double pr_iterate(Crystal *cr, const RefList *full) /* Calculate all gradients for this reflection */ for ( k=0; k<NUM_PARAMS; k++ ) { double gr; - gr = gradient(cr, k, refl); + gr = gradient(cr, k, refl, pmodel); gradients[k] = gr; } @@ -532,7 +524,7 @@ static double guide_dev(Crystal *cr, const RefList *full) } -void pr_refine(Crystal *cr, const RefList *full) +void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel) { double max_shift, dev; int i; @@ -557,9 +549,9 @@ void pr_refine(Crystal *cr, const RefList *full) cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - max_shift = pr_iterate(cr, full); + max_shift = pr_iterate(cr, full, pmodel); - update_partialities(cr); + update_partialities(cr, pmodel); if ( verbose ) { dev = guide_dev(cr, full); diff --git a/src/post-refinement.h b/src/post-refinement.h index fe171882..7b822938 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -40,6 +40,7 @@ #include "image.h" #include "utils.h" #include "crystal.h" +#include "geometry.h" /* Refineable parameters */ @@ -59,10 +60,11 @@ enum { }; -extern void pr_refine(Crystal *cr, const RefList *full); +extern void pr_refine(Crystal *cr, const RefList *full, PartialityModel pmodel); /* Exported so it can be poked by tests/pr_gradient_check */ -extern double gradient(Crystal *cr, int k, Reflection *refl); +extern double gradient(Crystal *cr, int k, Reflection *refl, + PartialityModel pmodel); #endif /* POST_REFINEMENT_H */ |