From 6815bab588d3e995bd1319e23614372aa1a07f0c Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 26 Apr 2011 15:13:42 +0200 Subject: Progress with post refinement... --- src/geometry.c | 29 +++++++++++++++++++ src/hrs-scaling.c | 15 ++++++---- src/partial_sim.c | 13 +++++++-- src/partialator.c | 41 +++----------------------- src/post-refinement.c | 79 +++++++++++++++++++++++++++++++++++++++++++++++---- 5 files changed, 126 insertions(+), 51 deletions(-) (limited to 'src') diff --git a/src/geometry.c b/src/geometry.c index b3d4c8ee..db261e6b 100644 --- a/src/geometry.c +++ b/src/geometry.c @@ -292,6 +292,34 @@ double integrate_all(struct image *image, RefList *reflections) } +/* Decide which reflections can be scaled */ +static int select_scalable_reflections(RefList *list) +{ + int n_scalable = 0; + + Reflection *refl; + RefListIterator *iter; + + for ( refl = first_refl(list, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + int scalable = 1; + double v; + + if ( get_partiality(refl) < 0.1 ) scalable = 0; + v = fabs(get_intensity(refl)); + if ( v < 0.1 ) scalable = 0; + + set_scalable(refl, scalable); + if ( scalable ) n_scalable++; + + } + + return n_scalable; +} + + /* Calculate partialities and apply them to the image's raw_reflections */ void update_partialities(struct image *image, const char *sym, int *n_expected, int *n_found, int *n_notfound) @@ -355,6 +383,7 @@ void update_partialities_and_asymm(struct image *image, const char *sym, * to get the list used for scaling and post refinement */ image->reflections = asymmetric_indices(image->raw_reflections, sym, obs); + select_scalable_reflections(image->reflections); /* Need these lists to work fast */ optimise_reflist(image->reflections); diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index 4b9b932d..cb0f8cde 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -338,12 +338,11 @@ static RefList *lsq_intensities(struct image *images, int n, for ( refl = find_refl(images[m].reflections, it->h, it->k, it->l); refl != NULL; - refl = next_found_refl(refl) ) { + refl = next_found_refl(refl) ) + { double p; - if ( !get_scalable(refl) ) continue; - p = get_partiality(refl); num += get_intensity(refl) * p * G; @@ -353,8 +352,14 @@ static RefList *lsq_intensities(struct image *images, int n, } - new = add_refl(full, it->h, it->k, it->l); - set_int(new, num/den); + if ( !isnan(num/den) ) { + new = add_refl(full, it->h, it->k, it->l); + set_int(new, num/den); + } else { + ERROR("Couldn't calculate LSQ full intensity for" + "%3i %3i %3i\n", it->h, it->k, it->l); + /* Doom is probably impending */ + } } diff --git a/src/partial_sim.c b/src/partial_sim.c index 824b1b5d..a9abb55b 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -39,10 +39,11 @@ static void mess_up_cell(UnitCell *cell) double cx, cy, cz; cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); - ax += 0.01 * ax; + ax += 0.005*ax; cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz); } + /* For each reflection in "partial", fill in what the intensity would be * according to "full" */ static void calculate_partials(RefList *partial, double osf, @@ -112,6 +113,7 @@ int main(int argc, char *argv[]) struct quaternion orientation; struct image image; FILE *ofh; + UnitCell *new; /* Long options */ const struct option longopts[] = { @@ -252,12 +254,17 @@ int main(int argc, char *argv[]) reflist_free(image.reflections); /* Alter the cell by a tiny amount */ - mess_up_cell(image.indexed_cell); image.filename = "(simulated 2)"; + new = rotate_cell(cell, deg2rad(0.001), deg2rad(0.0), 0.0); + cell_free(image.indexed_cell); + image.indexed_cell = new; - /* Write another chunk */ + /* Calculate new partials */ image.reflections = find_intersections(&image, image.indexed_cell, 0); calculate_partials(image.reflections, 0.5, full, sym); + + /* Give a slightly incorrect cell in the stream */ + //mess_up_cell(image.indexed_cell); write_chunk(ofh, &image, STREAM_INTEGRATED); reflist_free(image.reflections); diff --git a/src/partialator.c b/src/partialator.c index 54022e9e..9f931496 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -136,43 +136,12 @@ static void refine_all(struct image *images, int n_total_patterns, qargs.task_defaults = task_defaults; qargs.n = 0; qargs.n_done = 0; - qargs.n_total_patterns = n_total_patterns; - qargs.images = images; + /* FIXME: Not refining the first image, for now */ + qargs.n_total_patterns = n_total_patterns-1; + qargs.images = images+1; run_threads(nthreads, refine_image, get_image, done_image, - &qargs, n_total_patterns, 0, 0, 0); -} - - -/* Decide which reflections can be scaled */ -static void select_scalable_reflections(struct image *images, int n) -{ - int m; - int n_scalable = 0; - - for ( m=0; mreflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) { + + double G, p; + signed int h, k, l; + Reflection *full_version; + double I_full, I_partial; + + if ( !get_scalable(refl) ) continue; + + get_indices(refl, &h, &k, &l); + G = image->osf; + + full_version = find_refl(full, h, k, l); + assert(full_version != NULL); + + p = get_partiality(refl); + I_partial = get_intensity(refl); + I_full = get_intensity(full_version); + //STATUS("%3i %3i %3i %5.2f %5.2f %5.2f %5.2f\n", + // h, k, l, G, p, I_partial, I_full); + + dev += pow(I_partial - p*G*I_full, 2.0); + + } + + return dev; +} + + void pr_refine(struct image *image, const RefList *full, const char *sym) { double max_shift; int i; + double ax, ay, az; + double bx, by, bz; + double cx, cy, cz; + UnitCell *cell = image->indexed_cell; + double shval, origval; + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); + shval = 0.001*ax; + origval = ax; + + for ( i=-10; i<=10; i++ ) { + + double dev; + + cell_get_cartesian(cell, &ax, &ay, &az, &bx, + &by, &bz, &cx, &cy, &cz); + ax = origval + (double)i*shval; + cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz); + + update_partialities_and_asymm(image, sym, + NULL, NULL, NULL, NULL); + + dev = mean_partial_dev(image, full, sym); + STATUS("%i %e %e\n", i, ax, dev); + + } + return; i = 0; do { -- cgit v1.2.3