From c3c5ba136d1a4de7ab053f38724ece8c415b110d Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 26 Apr 2011 16:23:29 +0200 Subject: Work in progress on post refinement --- src/hrs-scaling.c | 3 ++- src/partial_sim.c | 6 ++--- src/post-refinement.c | 74 ++++++++++++++++++++++++++++++++++----------------- 3 files changed, 55 insertions(+), 28 deletions(-) diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c index cb0f8cde..fcf7efae 100644 --- a/src/hrs-scaling.c +++ b/src/hrs-scaling.c @@ -401,7 +401,8 @@ RefList *scale_intensities(struct image *images, int n, const char *sym, do { max_shift = iterate_scale(images, n, obs, sym, cref); - STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift); + STATUS("Scaling iteration %2i: max shift = %5.2f\n", + i, max_shift); i++; normalise_osfs(images, n); diff --git a/src/partial_sim.c b/src/partial_sim.c index a9abb55b..2c2493e9 100644 --- a/src/partial_sim.c +++ b/src/partial_sim.c @@ -39,7 +39,7 @@ 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.005*ax; + ax += 0.008*ax; cell_set_cartesian(cell, ax, ay, az, bx, by, bz, cx, cy, cz); } @@ -255,7 +255,7 @@ int main(int argc, char *argv[]) /* Alter the cell by a tiny amount */ image.filename = "(simulated 2)"; - new = rotate_cell(cell, deg2rad(0.001), deg2rad(0.0), 0.0); + new = rotate_cell(cell, deg2rad(1.0), deg2rad(0.0), 0.0); cell_free(image.indexed_cell); image.indexed_cell = new; @@ -264,7 +264,7 @@ int main(int argc, char *argv[]) calculate_partials(image.reflections, 0.5, full, sym); /* Give a slightly incorrect cell in the stream */ - //mess_up_cell(image.indexed_cell); + mess_up_cell(image.indexed_cell); write_chunk(ofh, &image, STREAM_INTEGRATED); reflist_free(image.reflections); diff --git a/src/post-refinement.c b/src/post-refinement.c index 7d851c17..20b8eeee 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -271,57 +271,63 @@ static double pr_iterate(struct image *image, const RefList *full, refl != NULL; refl = next_refl(refl, iter) ) { - signed int hind, kind, lind; signed int ha, ka, la; double I_full, delta_I; double I_partial; int k; double p; Reflection *match; - - get_indices(refl, &hind, &kind, &lind); + double gradients[NUM_PARAMS]; if ( !get_scalable(refl) ) continue; - /* Actual measurement of this reflection from this pattern? */ - I_partial = get_intensity(refl); - - get_asymm(hind, kind, lind, &ha, &ka, &la, sym); + /* Find the full version */ + get_indices(refl, &ha, &ka, &la); match = find_refl(full, ha, ka, la); - assert(match != NULL); + assert(match != NULL); /* Never happens because all scalable + * reflections had their LSQ intensities + * calculated in lsq_intensities(). */ I_full = get_intensity(match); + + /* Actual measurement of this reflection from this pattern? */ + I_partial = get_intensity(refl); p = get_partiality(refl); - delta_I = I_partial - (p * I_full / image->osf); + delta_I = I_partial - (p * image->osf * I_full); + + /* Calculate all gradients for this reflection */ + for ( k=0; kprofile_radius); + gradients[k] = gr; + } for ( k=0; kprofile_radius) - * gradient(image, k, refl, - image->profile_radius); + M_c = gradients[g] * gradients[k]; M_c *= pow(I_full, 2.0); + M_curr = gsl_matrix_get(M, g, k); gsl_matrix_set(M, g, k, M_curr + M_c); } gr = gradient(image, k, refl, image->profile_radius); v_c = delta_I * I_full * gr; - gsl_vector_set(v, k, v_c); + v_curr = gsl_vector_get(v, k); + gsl_vector_set(v, k, v_curr + v_c); } } - //show_matrix_eqn(M, v, NUM_PARAMS); + show_matrix_eqn(M, v, NUM_PARAMS); shifts = gsl_vector_alloc(NUM_PARAMS); gsl_linalg_HH_solve(M, v, shifts); @@ -381,15 +387,15 @@ static double mean_partial_dev(struct image *image, } -void pr_refine(struct image *image, const RefList *full, const char *sym) +static void plot_curve(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; + int i; cell_get_cartesian(cell, &ax, &ay, &az, &bx, &by, &bz, &cx, &cy, &cz); shval = 0.001*ax; @@ -411,12 +417,32 @@ void pr_refine(struct image *image, const RefList *full, const char *sym) STATUS("%i %e %e\n", i, ax, dev); } - return; +} + + +void pr_refine(struct image *image, const RefList *full, const char *sym) +{ + double max_shift; + int i; + + /* FIXME: This is for debugging */ + //plot_curve(image, full, sym); + //return; i = 0; do { + + double dev; + max_shift = pr_iterate(image, full, sym); - STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift); + update_partialities_and_asymm(image, sym, + NULL, NULL, NULL, NULL); + + dev = mean_partial_dev(image, full, sym); + STATUS("PR Iteration %2i: max shift = %5.2f dev = %5.2f\n", + i, max_shift, dev); + i++; + } while ( (max_shift > 0.01) && (i < MAX_CYCLES) ); } -- cgit v1.2.3