Commit efd6562f authored by Thomas White's avatar Thomas White
Browse files

Separate "refinable" and "scalable" concepts

parent b6742976
......@@ -23,6 +23,7 @@ get_symmetric_indices
get_intensity
get_partial
get_scalable
get_refinable
get_redundancy
get_sum_squared_dev
get_esd_intensity
......@@ -32,6 +33,7 @@ set_detector_pos
set_partial
set_int
set_scalable
set_refinable
set_redundancy
set_sum_squared_dev
set_esd_intensity
......
......@@ -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);
}
......@@ -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 */
......@@ -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);
......
......@@ -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++;
......
......@@ -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)
......@@ -381,6 +384,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)
......@@ -533,6 +548,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
......
......@@ -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);
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment