From 761da16096472d157b2dd9c871e81ab17cdef12d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Wed, 6 Feb 2013 18:59:40 +0100 Subject: Stuff for partialator --- src/hrs-scaling.c | 77 ++++++++++++++++++++++++++------------------------- src/partialator.c | 3 +- src/post-refinement.c | 64 +++++++++++++++++++++++------------------- src/post-refinement.h | 2 +- src/scaling-report.c | 70 +++++++++++++++++++++++++--------------------- 5 files changed, 116 insertions(+), 100 deletions(-) (limited to 'src') diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index f86ae0bb..0498a532 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -57,7 +57,7 @@ struct scale_queue_args { RefList *reference; - struct image *images; + Crystal **crystals; int n_started; double max_shift; }; @@ -65,7 +65,7 @@ struct scale_queue_args struct scale_worker_args { - struct image *image; + Crystal *crystal; double shift; RefList *reference; }; @@ -79,7 +79,7 @@ static void *create_scale_job(void *vqargs) wargs = malloc(sizeof(struct scale_worker_args)); wargs->reference = qargs->reference; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -88,20 +88,21 @@ static void *create_scale_job(void *vqargs) static void run_scale_job(void *vwargs, int cookie) { struct scale_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *reference = wargs->reference; Reflection *refl; RefListIterator *iter; double num = 0.0; double den = 0.0; double corr; + const double osf = crystal_get_osf(cr); - if ( image->pr_dud ) { + if ( crystal_get_user_flag(cr) ) { wargs->shift = 0.0; return; } - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -128,8 +129,8 @@ static void run_scale_job(void *vwargs, int cookie) Ihl = get_intensity(refl) / get_partiality(refl); esd = get_esd_intensity(refl) / get_partiality(refl); - num += Ih * (Ihl/image->osf) / pow(esd/image->osf, 2.0); - den += pow(Ih, 2.0)/pow(esd/image->osf, 2.0); + num += Ih * (Ihl/osf) / pow(esd/osf, 2.0); + den += pow(Ih, 2.0)/pow(esd/osf, 2.0); } @@ -138,7 +139,7 @@ static void run_scale_job(void *vwargs, int cookie) corr = num / den; if ( !isnan(corr) && !isinf(corr) ) { - image->osf *= corr; + crystal_set_osf(cr, osf*corr); } wargs->shift = fabs(corr-1.0); @@ -155,7 +156,7 @@ static void finalise_scale_job(void *vqargs, void *vwargs) } -static double iterate_scale(struct image *images, int n, RefList *reference, +static double iterate_scale(Crystal **crystals, int n, RefList *reference, int n_threads) { struct scale_queue_args qargs; @@ -164,7 +165,7 @@ static double iterate_scale(struct image *images, int n, RefList *reference, qargs.reference = reference; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; qargs.max_shift = 0.0; run_threads(n_threads, run_scale_job, create_scale_job, @@ -178,14 +179,14 @@ struct merge_queue_args { RefList *full; pthread_mutex_t full_lock; - struct image *images; + Crystal **crystals; int n_started; }; struct merge_worker_args { - struct image *image; + Crystal *crystal; RefList *full; pthread_mutex_t *full_lock; }; @@ -200,7 +201,7 @@ static void *create_merge_job(void *vqargs) wargs->full = qargs->full; wargs->full_lock = &qargs->full_lock; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -209,17 +210,17 @@ static void *create_merge_job(void *vqargs) static void run_merge_job(void *vwargs, int cookie) { struct merge_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *full = wargs->full; Reflection *refl; RefListIterator *iter; double G; - if ( image->pr_dud ) return; + if ( crystal_get_user_flag(cr)) return; - G = image->osf; + G = crystal_get_osf(cr); - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -272,7 +273,7 @@ static void finalise_merge_job(void *vqargs, void *vwargs) } -static RefList *lsq_intensities(struct image *images, int n, int n_threads) +static RefList *lsq_intensities(Crystal **crystals, int n, int n_threads) { RefList *full; struct merge_queue_args qargs; @@ -283,7 +284,7 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads) qargs.full = full; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; pthread_mutex_init(&qargs.full_lock, NULL); run_threads(n_threads, run_merge_job, create_merge_job, @@ -308,14 +309,14 @@ static RefList *lsq_intensities(struct image *images, int n, int n_threads) struct esd_queue_args { RefList *full; - struct image *images; + Crystal **crystals; int n_started; }; struct esd_worker_args { - struct image *image; + Crystal *crystal; RefList *full; }; @@ -328,7 +329,7 @@ static void *create_esd_job(void *vqargs) wargs = malloc(sizeof(struct esd_worker_args)); wargs->full = qargs->full; - wargs->image = &qargs->images[qargs->n_started++]; + wargs->crystal = qargs->crystals[qargs->n_started++]; return wargs; } @@ -337,17 +338,17 @@ static void *create_esd_job(void *vqargs) static void run_esd_job(void *vwargs, int cookie) { struct esd_worker_args *wargs = vwargs; - struct image *image = wargs->image; + Crystal *cr = wargs->crystal; RefList *full = wargs->full; Reflection *refl; RefListIterator *iter; double G; - if ( image->pr_dud ) return; + if ( crystal_get_user_flag(cr) ) return; - G = image->osf; + G = crystal_get_osf(cr); - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -383,7 +384,7 @@ static void finalise_esd_job(void *vqargs, void *vwargs) } -static void calculate_esds(struct image *images, int n, RefList *full, +static void calculate_esds(Crystal **crystals, int n, RefList *full, int n_threads, int min_red) { struct esd_queue_args qargs; @@ -392,7 +393,7 @@ static void calculate_esds(struct image *images, int n, RefList *full, qargs.full = full; qargs.n_started = 0; - qargs.images = images; + qargs.crystals = crystals; for ( refl = first_refl(full, &iter); refl != NULL; @@ -423,7 +424,7 @@ static void calculate_esds(struct image *images, int n, RefList *full, /* Scale the stack of images */ -RefList *scale_intensities(struct image *images, int n, RefList *gref, +RefList *scale_intensities(Crystal **crystals, int n, RefList *gref, int n_threads, int noscale) { int i; @@ -431,17 +432,17 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, RefList *full = NULL; const int min_redundancy = 3; - for ( i=0; i create an initial list to refine against */ if ( gref == NULL ) { - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } /* Iterate */ @@ -457,14 +458,14 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, reference = full; } - max_corr = iterate_scale(images, n, reference, n_threads); + max_corr = iterate_scale(crystals, n, reference, n_threads); //STATUS("Scaling iteration %2i: max correction = %5.2f\n", // i+1, max_corr); /* No reference -> generate list for next iteration */ if ( gref == NULL ) { reflist_free(full); - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } //show_scale_factors(images, n); @@ -474,10 +475,10 @@ RefList *scale_intensities(struct image *images, int n, RefList *gref, } while ( (max_corr > 0.01) && (i < MAX_CYCLES) ); if ( gref != NULL ) { - full = lsq_intensities(images, n, n_threads); + full = lsq_intensities(crystals, n, n_threads); } /* else we already did it */ - calculate_esds(images, n, full, n_threads, min_redundancy); + calculate_esds(crystals, n, full, n_threads, min_redundancy); return full; } diff --git a/src/partialator.c b/src/partialator.c index 35e25320..30f9bfe4 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -506,6 +506,7 @@ int main(int argc, char *argv[]) crystals = crystals_new; crystals[n_crystals] = crystal_new(); cr = crystals[n_crystals]; + crystal_set_image(cr, cur); /* Fill in initial estimates of stuff */ crystal_set_osf(cr, 1.0); @@ -516,7 +517,7 @@ int main(int argc, char *argv[]) as = asymmetric_indices(crystal_get_reflections(cr), sym); crystal_set_reflections(cr, as); - update_partialities(cur, cr); + update_partialities(cr); nobs += select_scalable_reflections(as, reference); diff --git a/src/post-refinement.c b/src/post-refinement.c index 7410d931..31558ec7 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -3,11 +3,11 @@ * * Post refinement * - * Copyright © 2012 Deutsches Elektronen-Synchrotron DESY, - * a research centre of the Helmholtz Association. + * Copyright © 2012-2013 Deutsches Elektronen-Synchrotron DESY, + * a research centre of the Helmholtz Association. * * Authors: - * 2010-2012 Thomas White + * 2010-2013 Thomas White * * This file is part of CrystFEL. * @@ -88,7 +88,7 @@ static double partiality_rgradient(double r, double profile_radius) /* Return the gradient of parameter 'k' given the current status of 'image'. */ -double gradient(struct image *image, int k, Reflection *refl, double r) +double gradient(Crystal *cr, int k, Reflection *refl, double r) { double ds, azix, aziy; double ttlow, tthigh, tt; @@ -103,17 +103,18 @@ double gradient(struct image *image, int k, Reflection *refl, double r) int clamp_low, clamp_high; double klow, khigh; double gr; + struct image *image = crystal_get_image(cr); get_symmetric_indices(refl, &hs, &ks, &ls); - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, - &bsx, &bsy, &bsz, - &csx, &csy, &csz); + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, + &bsx, &bsy, &bsz, + &csx, &csy, &csz); xl = hs*asx + ks*bsx + ls*csx; yl = hs*asy + ks*bsy + ls*csy; zl = hs*asz + ks*bsz + ls*csz; - ds = 2.0 * resolution(image->indexed_cell, hs, ks, ls); + ds = 2.0 * resolution(crystal_get_cell(cr), hs, ks, ls); get_partial(refl, &r1, &r2, &p, &clamp_low, &clamp_high); klow = 1.0/(image->lambda - image->lambda*image->bw/2.0); @@ -237,8 +238,11 @@ static void apply_cell_shift(UnitCell *cell, int k, double shift) /* Apply the given shift to the 'k'th parameter of 'image'. */ -static void apply_shift(struct image *image, int k, double shift) +static void apply_shift(Crystal *cr, int k, double shift) { + double t; + struct image *image = crystal_get_image(cr); + switch ( k ) { case REF_DIV : @@ -251,7 +255,9 @@ static void apply_shift(struct image *image, int k, double shift) break; case REF_R : - image->profile_radius += shift; + t = crystal_get_profile_radius(cr); + t += shift; + crystal_set_profile_radius(cr, t); break; case REF_ASX : @@ -263,7 +269,7 @@ static void apply_shift(struct image *image, int k, double shift) case REF_CSX : case REF_CSY : case REF_CSZ : - apply_cell_shift(image->indexed_cell, k, shift); + apply_cell_shift(crystal_get_cell(cr), k, shift); break; default : @@ -357,7 +363,7 @@ static gsl_vector *solve_svd(gsl_vector *v, gsl_matrix *M) /* Perform one cycle of post refinement on 'image' against 'full' */ -static double pr_iterate(struct image *image, const RefList *full) +static double pr_iterate(Crystal *cr, const RefList *full) { gsl_matrix *M; gsl_vector *v; @@ -369,7 +375,7 @@ static double pr_iterate(struct image *image, const RefList *full) double max_shift; int nref = 0; - reflections = image->reflections; + reflections = crystal_get_reflections(cr); M = gsl_matrix_calloc(NUM_PARAMS, NUM_PARAMS); v = gsl_vector_calloc(NUM_PARAMS); @@ -387,6 +393,7 @@ static double pr_iterate(struct image *image, const RefList *full) double p; Reflection *match; double gradients[NUM_PARAMS]; + const double osf = crystal_get_osf(cr); if ( !get_refinable(refl) ) continue; @@ -398,7 +405,7 @@ static double pr_iterate(struct image *image, const RefList *full) " is it marked as refinable?\n", ha, ka, la); continue; } - I_full = image->osf * get_intensity(match); + I_full = osf * get_intensity(match); /* Actual measurement of this reflection from this pattern? */ I_partial = get_intensity(refl); @@ -412,7 +419,8 @@ static double pr_iterate(struct image *image, const RefList *full) /* Calculate all gradients for this reflection */ for ( k=0; kprofile_radius); + gr = gradient(cr, k, refl, + crystal_get_profile_radius(cr)); gradients[k] = gr; } @@ -429,7 +437,7 @@ static double pr_iterate(struct image *image, const RefList *full) if ( g > k ) continue; M_c = gradients[g] * gradients[k]; - M_c *= w * pow(image->osf * I_full, 2.0); + M_c *= w * pow(osf * I_full, 2.0); M_curr = gsl_matrix_get(M, k, g); gsl_matrix_set(M, k, g, M_curr + M_c); @@ -450,7 +458,7 @@ static double pr_iterate(struct image *image, const RefList *full) //STATUS("%i reflections went into the equations.\n", nref); if ( nref == 0 ) { - image->pr_dud = 1; + crystal_set_user_flag(cr, 1); gsl_matrix_free(M); gsl_vector_free(v); return 0.0; @@ -462,7 +470,7 @@ static double pr_iterate(struct image *image, const RefList *full) for ( param=0; param max_shift ) max_shift = fabs(shift); } @@ -480,7 +488,7 @@ static double pr_iterate(struct image *image, const RefList *full) } -static double guide_dev(struct image *image, const RefList *full) +static double guide_dev(Crystal *cr, const RefList *full) { double dev = 0.0; @@ -488,7 +496,7 @@ static double guide_dev(struct image *image, const RefList *full) Reflection *refl; RefListIterator *iter; - for ( refl = first_refl(image->reflections, &iter); + for ( refl = first_refl(crystal_get_reflections(cr), &iter); refl != NULL; refl = next_refl(refl, iter) ) { @@ -508,7 +516,7 @@ static double guide_dev(struct image *image, const RefList *full) * scale_intensities() might not yet have been called, so the * full version may not have been calculated yet. */ - G = image->osf; + G = crystal_get_osf(cr); p = get_partiality(refl); I_partial = get_intensity(refl); I_full = get_intensity(full_version); @@ -524,21 +532,21 @@ static double guide_dev(struct image *image, const RefList *full) } -void pr_refine(struct image *image, const RefList *full) +void pr_refine(Crystal *cr, const RefList *full) { double max_shift, dev; int i; const int verbose = 0; if ( verbose ) { - dev = guide_dev(image, full); + dev = guide_dev(cr, full); STATUS("\n"); /* Deal with progress bar */ STATUS("Before iteration: dev = %10.5e\n", dev); } i = 0; - image->pr_dud = 0; + crystal_set_user_flag(cr, 0); do { double asx, asy, asz; @@ -546,15 +554,15 @@ void pr_refine(struct image *image, const RefList *full) double csx, csy, csz; double dev; - cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, + cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - max_shift = pr_iterate(image, full); + max_shift = pr_iterate(cr, full); - update_partialities(image); + update_partialities(cr); if ( verbose ) { - dev = guide_dev(image, full); + dev = guide_dev(cr, full); STATUS("PR Iteration %2i: max shift = %10.2f" " dev = %10.5e\n", i+1, max_shift, dev); } diff --git a/src/post-refinement.h b/src/post-refinement.h index 637c0bb0..52e70b2b 100644 --- a/src/post-refinement.h +++ b/src/post-refinement.h @@ -62,7 +62,7 @@ enum { extern void pr_refine(Crystal *cr, const RefList *full); /* Exported so it can be poked by tests/pr_gradient_check */ -extern double gradient(struct image *image, int k, Reflection *refl, double r); +extern double gradient(Crystal *cr, int k, Reflection *refl, double r); #endif /* POST_REFINEMENT_H */ diff --git a/src/scaling-report.c b/src/scaling-report.c index ba425b96..3de82517 100644 --- a/src/scaling-report.c +++ b/src/scaling-report.c @@ -183,7 +183,7 @@ static void plot_point(cairo_t *cr, double g_width, double g_height, } -static void partiality_graph(cairo_t *cr, const struct image *images, int n, +static void partiality_graph(cairo_t *cr, Crystal **crystals, int n, RefList *full) { const double g_width = 200.0; @@ -219,7 +219,7 @@ static void partiality_graph(cairo_t *cr, const struct image *images, int n, num_nondud = 0; for ( i=0; i osf_max ) osf_max = osf; } osf_max = ceil(osf_max+osf_max/10000.0); @@ -473,9 +477,9 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images, for ( i=0; i= osf_low[b]) @@ -544,8 +548,8 @@ static void scale_factor_histogram(cairo_t *cr, const struct image *images, } -static void intensity_histogram(cairo_t *cr, const struct image *images, - int n, RefList *full, +static void intensity_histogram(cairo_t *cr, Crystal **crystals, int n, + RefList *full, signed int h, signed int k, signed int l) { int f_max; @@ -582,12 +586,13 @@ static void intensity_histogram(cairo_t *cr, const struct image *images, for ( i=0; icr); cairo_translate(sr->cr, 480.0, 350.0); - scale_factor_histogram(sr->cr, images, n, + scale_factor_histogram(sr->cr, crystals, n, "Distribution of overall scale factors"); cairo_restore(sr->cr); @@ -830,7 +836,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_save(sr->cr); cairo_translate(sr->cr, 70.0, 330.0); - partiality_graph(sr->cr, images, n, full); + partiality_graph(sr->cr, crystals, n, full); cairo_save(sr->cr); cairo_move_to(sr->cr, 0.0, 0.0); @@ -841,7 +847,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_stroke(sr->cr); cairo_set_dash(sr->cr, NULL, 0, 0.0); cairo_translate(sr->cr, 0.0, -150.0); - partiality_histogram(sr->cr, images, n, full, 1, 0); + partiality_histogram(sr->cr, crystals, n, full, 1, 0); cairo_restore(sr->cr); cairo_save(sr->cr); @@ -854,7 +860,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_set_dash(sr->cr, NULL, 0, 0.0); cairo_translate(sr->cr, 230.0, 200.0); cairo_rotate(sr->cr, -M_PI_2); - partiality_histogram(sr->cr, images, n, full, 0, 1); + partiality_histogram(sr->cr, crystals, n, full, 0, 1); cairo_restore(sr->cr); cairo_restore(sr->cr); @@ -873,7 +879,7 @@ void sr_iteration(SRContext *sr, int iteration, struct image *images, int n, cairo_save(sr->cr); cairo_translate(sr->cr, 400.0+140.0*x, 60.0+80.0*y); - intensity_histogram(sr->cr, images, n, full, + intensity_histogram(sr->cr, crystals, n, full, sr->ms_h[i], sr->ms_k[i], sr->ms_l[i]); cairo_restore(sr->cr); -- cgit v1.2.3