diff options
Diffstat (limited to 'src')
| -rw-r--r-- | src/geometry.c | 48 | ||||
| -rw-r--r-- | src/geometry.h | 8 | ||||
| -rw-r--r-- | src/partialator.c | 45 | ||||
| -rw-r--r-- | src/post-refinement.c | 49 | ||||
| -rw-r--r-- | src/reflist.c | 35 | ||||
| -rw-r--r-- | src/reflist.h | 2 |
6 files changed, 131 insertions, 56 deletions
diff --git a/src/geometry.c b/src/geometry.c index e2454adf..7a11b110 100644 --- a/src/geometry.c +++ b/src/geometry.c @@ -275,11 +275,10 @@ RefList *find_intersections(struct image *image, UnitCell *cell) } -/* Calculate partialities and apply them to the image's raw_reflections, - * while adding to a ReflItemList of the currentl scalable (asymmetric) - * reflections. */ -void update_partialities(struct image *image, const char *sym, - int *n_expected, int *n_found, int *n_notfound) +/* Predict reflections in "image" */ +void predict_corresponding_reflections(struct image *image, const char *sym, + int *n_expected, int *n_found, + int *n_notfound) { Reflection *refl; RefListIterator *iter; @@ -346,3 +345,42 @@ void update_partialities(struct image *image, const char *sym, reflist_free(predicted); } + + +/* Calculate partialities and apply them to the image's raw_reflections */ +void update_partialities(struct image *image, const char *sym) +{ + Reflection *refl; + RefListIterator *iter; + RefList *predicted; + + predicted = find_intersections(image, image->indexed_cell); + + for ( refl = first_refl(image->reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + + Reflection *p_peak; + double r1, r2, p; + signed int h, k, l; + int clamp1, clamp2; + + /* Get predicted indices and location */ + get_symmetric_indices(refl, &h, &k, &l); + + /* Look for this reflection in the pattern */ + p_peak = find_refl(predicted, h, k, l); + if ( p_peak == NULL ) { + set_partial(refl, 0.0, 0.0, 0.0, -1, +1); + continue; + } else { + /* Transfer partiality stuff */ + get_partial(p_peak, &r1, &r2, &p, &clamp1, &clamp2); + set_partial(refl, r1, r2, p, clamp1, clamp2); + } + + } + + reflist_free(predicted); +} diff --git a/src/geometry.h b/src/geometry.h index 892ff747..efda6c87 100644 --- a/src/geometry.h +++ b/src/geometry.h @@ -21,7 +21,11 @@ extern RefList *find_intersections(struct image *image, UnitCell *cell); -extern void update_partialities(struct image *image, const char *sym, - int *n_expected, int *n_found, int *n_notfound); +extern void predict_corresponding_reflections(struct image *image, + const char *sym, int *n_expected, + int *n_found, int *n_notfound); + +extern void update_partialities(struct image *image, const char *sym); + #endif /* GEOMETRY_H */ diff --git a/src/partialator.c b/src/partialator.c index 167f736d..43c4bda0 100644 --- a/src/partialator.c +++ b/src/partialator.c @@ -191,6 +191,35 @@ static int select_scalable_reflections(RefList *list, ReflItemList *sc_l) } +static void select_reflections_for_refinement(struct image *images, int n, + ReflItemList *scalable) +{ + int i; + + for ( i=0; i<n; i++ ) { + + Reflection *refl; + RefListIterator *iter; + + for ( refl = first_refl(images[i].reflections, &iter); + refl != NULL; + refl = next_refl(refl, iter) ) + { + signed int h, k, l; + int sc; + + get_indices(refl, &h, &k, &l); + sc = get_scalable(refl); + + if ( sc && find_item(scalable, h, k, l) ) { + set_refinable(refl, 1); + } + } + + } +} + + int main(int argc, char *argv[]) { int c; @@ -397,8 +426,8 @@ int main(int argc, char *argv[]) reflist_free(cur->reflections); cur->reflections = as; - update_partialities(cur, sym, - &n_expected, &n_found, &n_notfound); + predict_corresponding_reflections(cur, sym, &n_expected, + &n_found, &n_notfound); nobs += select_scalable_reflections(cur->reflections, scalable); @@ -419,6 +448,8 @@ int main(int argc, char *argv[]) full = scale_intensities(images, n_usable_patterns, sym, scalable, cref, reference); + select_reflections_for_refinement(images, n_usable_patterns, scalable); + for ( i=0; i<num_items(scalable); i++ ) { Reflection *f; struct refl_item *it = get_item(scalable, i); @@ -459,6 +490,7 @@ int main(int argc, char *argv[]) FILE *fhg; FILE *fhp; char filename[1024]; + int j; STATUS("Post refinement cycle %i of %i\n", i+1, n_iter); @@ -484,9 +516,14 @@ int main(int argc, char *argv[]) nobs = 0; clear_items(scalable); - for ( i=0; i<n_usable_patterns; i++ ) { + for ( j=0; j<n_usable_patterns; j++ ) { + + struct image *cur = &images[j]; + + predict_corresponding_reflections(cur, sym, &n_expected, + &n_found, + &n_notfound); - struct image *cur = &images[i]; nobs += select_scalable_reflections(cur->reflections, scalable); diff --git a/src/post-refinement.c b/src/post-refinement.c index 7eef5dfa..aeafcddc 100644 --- a/src/post-refinement.c +++ b/src/post-refinement.c @@ -31,7 +31,7 @@ /* Maximum number of iterations of NLSq to do for each image per macrocycle. */ -#define MAX_CYCLES (5) +#define MAX_CYCLES (10) /* Returns dp/dr at "r" */ @@ -347,7 +347,7 @@ static double pr_iterate(struct image *image, const RefList *full, Reflection *match; double gradients[NUM_PARAMS]; - if ( !get_scalable(refl) ) continue; + if ( !get_refinable(refl) ) continue; /* Find the full version */ get_indices(refl, &ha, &ka, &la); @@ -484,66 +484,33 @@ void pr_refine(struct image *image, const RefList *full, const char *sym) double max_shift, dev; int i; const int verbose = 1; - int nexp, nfound, nnotfound; - - update_partialities(image, sym, &nexp, &nfound, &nnotfound); if ( verbose ) { dev = mean_partial_dev(image, full, sym); - STATUS("PR starting dev = %5.2f (%i out of %i found)\n", - dev, nfound, nexp); - } - - if ( (double)nfound/(double)nexp < 0.5 ) { - ERROR("Refusing to refine this image: %i out of %i found\n", - nfound, nexp); - return; + STATUS("PR starting dev = %5.2f\n", dev); } i = 0; + image->pr_dud = 0; do { double asx, asy, asz; double bsx, bsy, bsz; double csx, csy, csz; double dev; - int old_nexp, old_nfound; cell_get_reciprocal(image->indexed_cell, &asx, &asy, &asz, &bsx, &bsy, &bsz, &csx, &csy, &csz); - old_nexp = nexp; - old_nfound = nfound; max_shift = pr_iterate(image, full, sym); - update_partialities(image, sym, &nexp, &nfound, &nnotfound); + update_partialities(image, sym); if ( verbose ) { dev = mean_partial_dev(image, full, sym); - STATUS("PR Iteration %2i: max shift = %5.2f" - " dev = %5.2f (%i out of %i found)\n", - i+1, max_shift, dev, nfound, nexp); - } - - if ( (double)nfound / (double)nexp < 0.5 ) { - - if ( verbose ) { - ERROR("Bad refinement step - backtracking.\n"); - ERROR("I'll come back to this image later.\n"); - } - - cell_set_reciprocal(image->indexed_cell, asx, asy, asz, - bsx, bsy, bsz, csx, csy, csz); - - update_partialities(image, sym, - &nexp, &nfound, &nnotfound); - - image->pr_dud = 1; - - return; - - } else { - image->pr_dud = 0; + STATUS("PR Iteration %2i: max shift = %10.2f" + " dev = %10.5e\n", + i+1, max_shift, dev); } i++; diff --git a/src/reflist.c b/src/reflist.c index a8bfb4f6..19d5b1a4 100644 --- a/src/reflist.c +++ b/src/reflist.c @@ -67,6 +67,10 @@ struct _refldata { /* Non-zero if this reflection can be used for scaling */ int scalable; + /* Non-zero if this reflection should be used as a "guide star" for + * post refinement */ + int refinable; + /* Intensity */ double intensity; double esd_i; @@ -371,8 +375,7 @@ void get_partial(const Reflection *refl, double *r1, double *r2, double *p, * get_scalable: * @refl: A %Reflection * - * Returns: non-zero if this reflection was marked as useful for scaling and - * post refinement. + * Returns: non-zero if this reflection can be scaled. * **/ int get_scalable(const Reflection *refl) @@ -382,6 +385,19 @@ int get_scalable(const Reflection *refl) /** + * get_refinable: + * @refl: A %Reflection + * + * Returns: non-zero if this reflection can be used for post refinement. + * + **/ +int get_refinable(const Reflection *refl) +{ + return refl->data.refinable; +} + + +/** * get_redundancy: * @refl: A %Reflection * @@ -523,8 +539,7 @@ void set_int(Reflection *refl, double intensity) /** * set_scalable: * @refl: A %Reflection - * @scalable: Non-zero if this reflection was marked as useful for scaling and - * post refinement. + * @scalable: Non-zero if this reflection should be scaled. * **/ void set_scalable(Reflection *refl, int scalable) @@ -534,6 +549,18 @@ void set_scalable(Reflection *refl, int scalable) /** + * set_refinable: + * @refl: A %Reflection + * @scalable: Non-zero if this reflection can be used for post refinement. + * + **/ +void set_refinable(Reflection *refl, int refinable) +{ + refl->data.refinable = refinable; +} + + +/** * set_redundancy: * @refl: A %Reflection * @red: New redundancy for the reflection diff --git a/src/reflist.h b/src/reflist.h index ffdef55a..0d836edf 100644 --- a/src/reflist.h +++ b/src/reflist.h @@ -68,6 +68,7 @@ extern double get_intensity(const Reflection *refl); extern void get_partial(const Reflection *refl, double *r1, double *r2, double *p, int *clamp_low, int *clamp_high); extern int get_scalable(const Reflection *refl); +extern int get_refinable(const Reflection *refl); extern int get_redundancy(const Reflection *refl); extern double get_sum_squared_dev(const Reflection *refl); extern double get_esd_intensity(const Reflection *refl); @@ -81,6 +82,7 @@ extern void set_partial(Reflection *refl, double r1, double r2, double p, double clamp_low, double clamp_high); extern void set_int(Reflection *refl, double intensity); extern void set_scalable(Reflection *refl, int scalable); +extern void set_refinable(Reflection *refl, int refinable); extern void set_redundancy(Reflection *refl, int red); extern void set_sum_squared_dev(Reflection *refl, double dev); extern void set_esd_intensity(Reflection *refl, double esd); |