aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-07-04 17:50:56 +0200
committerThomas White <taw@physics.org>2012-02-22 15:27:31 +0100
commitefd6562f9156ddff2fe073c97b2ddbf25c45688e (patch)
tree1d0320634e5cacf9130900b6a013010503a34f41 /src
parentb67429762f02d906fdc3ab14da4577c958937679 (diff)
Separate "refinable" and "scalable" concepts
Diffstat (limited to 'src')
-rw-r--r--src/geometry.c48
-rw-r--r--src/geometry.h8
-rw-r--r--src/partialator.c45
-rw-r--r--src/post-refinement.c49
-rw-r--r--src/reflist.c35
-rw-r--r--src/reflist.h2
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);