From cef0a71eb385773fa2290dfc99de225948fef06a Mon Sep 17 00:00:00 2001 From: Thomas White Date: Fri, 15 Mar 2013 15:53:39 +0100 Subject: partialator: Add --model= --- libcrystfel/src/geometry.c | 23 +++++++++++ src/hrs-scaling.c | 101 ++++++++++++++++++++++++++++++++++----------- src/hrs-scaling.h | 9 ++-- src/partialator.c | 17 +++++--- 4 files changed, 116 insertions(+), 34 deletions(-) diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c index f457fc4f..cdc62ab5 100644 --- a/libcrystfel/src/geometry.c +++ b/libcrystfel/src/geometry.c @@ -362,6 +362,21 @@ RefList *select_intersections(struct image *image, Crystal *cryst) } +static void set_unity_partialities(Crystal *cryst) +{ + Reflection *refl; + RefListIterator *iter; + + for ( refl = first_refl(crystal_get_reflections(cryst), &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + set_partiality(refl, 1.0); + set_lorentz(refl, 1.0); + } +} + + /* Calculate partialities and apply them to the image's reflections */ void update_partialities(Crystal *cryst, PartialityModel pmodel) { @@ -373,6 +388,14 @@ void update_partialities(Crystal *cryst, PartialityModel pmodel) double csx, csy, csz; struct image *image = crystal_get_image(cryst); + if ( pmodel == PMODEL_UNITY ) { + /* It isn't strictly necessary to set the partialities to 1, + * because the scaling stuff will just a correction factor of + * 1 anyway. This is just to help things make sense. */ + set_unity_partialities(cryst); + return; + } + cell_get_reciprocal(crystal_get_cell(cryst), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index f60dfe25..3a51e6d9 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -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++]; @@ -109,7 +112,7 @@ static void run_scale_job(void *vwargs, int cookie) signed int h, k, l; double Ih, Ihl, esd; Reflection *r; - double L, p; + double corr; if ( !get_scalable(refl) ) continue; @@ -126,11 +129,21 @@ static void run_scale_job(void *vwargs, int cookie) Ih = get_intensity(r); } - /* If you change this, be sure to change run_merge_job() too */ - p = get_partiality(refl); - L = get_lorentz(refl); - Ihl = get_intensity(refl) / (L*p); - esd = get_esd_intensity(refl) / (L*p); + /* If you change this, be sure to also change + * run_merge_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; + + } + Ihl = get_intensity(refl) / corr; + esd = get_esd_intensity(refl) / corr; num += Ih * (Ihl/osf) / pow(esd/osf, 2.0); den += pow(Ih, 2.0)/pow(esd/osf, 2.0); @@ -160,7 +173,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; @@ -170,6 +183,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); @@ -184,6 +198,7 @@ struct merge_queue_args pthread_mutex_t full_lock; Crystal **crystals; int n_started; + PartialityModel pmodel; }; @@ -192,6 +207,7 @@ struct merge_worker_args Crystal *crystal; RefList *full; pthread_mutex_t *full_lock; + PartialityModel pmodel; }; @@ -203,6 +219,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++]; @@ -231,7 +248,7 @@ static void run_merge_job(void *vwargs, int cookie) signed int h, k, l; double num, den; int red; - double Ihl, esd, L, p; + double Ihl, esd, corr; if ( !get_scalable(refl) ) continue; @@ -254,11 +271,22 @@ static void run_merge_job(void *vwargs, int cookie) red = get_redundancy(f); } - /* If you change this, be sure to change run_scale_job() too */ - p = get_partiality(refl); - L = get_lorentz(refl); - Ihl = get_intensity(refl) / (L*p); - esd = get_esd_intensity(refl) / (L*p); + /* 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;; + + } + + 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); @@ -278,7 +306,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; @@ -290,6 +319,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, @@ -316,6 +346,7 @@ struct esd_queue_args RefList *full; Crystal **crystals; int n_started; + PartialityModel pmodel; }; @@ -323,6 +354,7 @@ struct esd_worker_args { Crystal *crystal; RefList *full; + PartialityModel pmodel; }; @@ -333,6 +365,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++]; @@ -360,7 +393,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; @@ -372,8 +405,23 @@ 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;; + + } + + Ih = get_intensity(f); - Ihl = get_intensity(refl) / (get_partiality(refl) * G); + Ihl = get_intensity(refl) / (corr * G); num += pow(Ihl - Ih, 2.0); @@ -390,7 +438,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; @@ -399,6 +447,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; @@ -430,7 +479,7 @@ 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 i; double max_corr; @@ -440,14 +489,15 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, for ( i=0; i 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 */ @@ -463,14 +513,15 @@ 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); @@ -480,10 +531,10 @@ RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, } 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..8ae12380 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 + * 2010-2013 Thomas White * * This file is part of CrystFEL. * @@ -37,10 +37,11 @@ #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); #endif /* HRS_SCALING_H */ diff --git a/src/partialator.c b/src/partialator.c index 079e078e..bdacd62b 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -79,6 +79,7 @@ static void show_help(const char *s) " -r, --reference= Refine images against reflections in ,\n" " instead of taking the mean of the intensity\n" " estimates.\n" +" -m, --model= Specify partiality model.\n" "\n" " -j Run analyses in parallel.\n"); } @@ -140,13 +141,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; @@ -329,7 +335,8 @@ int main(int argc, char *argv[]) {"iterations", 1, NULL, 'n'}, {"no-scale", 0, &noscale, 1}, {"reference", 1, NULL, 'r'}, - {"partiality", 1, NULL, 'm'}, + {"model", 1, NULL, 'm'}, + {0, 0, NULL, 0} }; @@ -574,7 +581,7 @@ int main(int argc, char *argv[]) STATUS("\nPerforming initial scaling.\n"); if ( noscale ) STATUS("Scale factors fixed at 1.\n"); full = scale_intensities(crystals, n_crystals, reference, - nthreads, noscale); + nthreads, noscale, pmodel); sr = sr_titlepage(crystals, n_crystals, "scaling-report.pdf", infile, cmdline); @@ -597,7 +604,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