Skip to content
GitLab
Commit 7a7cf035 authored by Thomas White's avatar Thomas White
Browse files

Fix the choice of "guiding" reflections for PR

parent d5337c0c
Hide whitespace changes
Inline Side-by-side
src/geometry.c
View file @ 7a7cf035
......@@ -347,7 +347,7 @@ void predict_corresponding_reflections(struct image *image, const char *sym,
/* Calculate partialities and apply them to the image's raw_reflections */
void update_partialities(struct image *image, const char *sym)
void update_partialities(struct image *image)
{
Reflection *refl;
RefListIterator *iter;
......
src/geometry.h
View file @ 7a7cf035
......@@ -25,7 +25,7 @@ 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);
extern void update_partialities(struct image *image);
#endif /* GEOMETRY_H */
src/hrs-scaling.c
View file @ 7a7cf035
......@@ -35,7 +35,7 @@
#define MAX_CYCLES (50)
char *find_common_reflections(struct image *images, int n)
static char *find_common_reflections(struct image *images, int n)
{
int i;
char *cref;
......@@ -54,6 +54,8 @@ char *find_common_reflections(struct image *images, int n)
signed int h, k, l;
int j;
if ( !get_scalable(refl) ) continue;
get_indices(refl, &h, &k, &l);
for ( j=0; j<i; j++ ) {
......@@ -84,8 +86,8 @@ static void s_uhavha(signed int hat, signed int kat, signed int lat,
for ( refl = find_refl(reflections, hat, kat, lat);
refl != NULL;
refl = next_found_refl(refl) ) {
refl = next_found_refl(refl) )
{
double ic, sigi;
if ( !get_scalable(refl) ) continue;
......@@ -95,7 +97,6 @@ static void s_uhavha(signed int hat, signed int kat, signed int lat,
uha_val += 1.0 / pow(sigi, 2.0);
vha_val += ic / pow(sigi, 2.0);
}
*uha = uha_val;
......@@ -104,7 +105,7 @@ static void s_uhavha(signed int hat, signed int kat, signed int lat,
static void s_uhvh(struct image *images, int n,
signed int h, signed int k, signed int l, const char *sym,
signed int h, signed int k, signed int l,
double *uhp, double *vhp)
{
int a;
......@@ -212,51 +213,52 @@ static gsl_vector *solve_diagonal(gsl_vector *v, gsl_matrix *M)
}
static double iterate_scale(struct image *images, int n,
ReflItemList *obs, const char *sym, char *cref,
double *reference)
static double iterate_scale(struct image *images, int n, RefList *scalable,
char *cref, RefList *reference)
{
gsl_matrix *M;
gsl_vector *v;
gsl_vector *shifts;
double max_shift;
int n_ref;
double *uh_arr;
double *vh_arr;
double *uha_arr;
double *vha_arr;
int h; /* Reflection index */
int frame;
int refidx;
Reflection *refl;
RefListIterator *iter;
M = gsl_matrix_calloc(n, n);
v = gsl_vector_calloc(n);
n_ref = num_items(obs);
uh_arr = new_list_intensity();
vh_arr = new_list_intensity();
for ( h=0; h<n_ref; h++ ) {
for ( refl = first_refl(scalable, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
double uh, vh;
struct refl_item *it = get_item(obs, h);
signed int h, k, l;
s_uhvh(images, n, it->h, it->k, it->l, sym, &uh, &vh);
get_indices(refl, &h, &k, &l);
set_intensity(uh_arr, it->h, it->k, it->l, uh);
set_intensity(vh_arr, it->h, it->k, it->l, vh);
s_uhvh(images, n, h, k, l, &uh, &vh);
set_intensity(uh_arr, h, k, l, uh);
set_intensity(vh_arr, h, k, l, vh);
}
uha_arr = malloc(n*sizeof(double));
vha_arr = malloc(n*sizeof(double));
for ( refidx=0; refidx<n_ref; refidx++ ) {
for ( refl = first_refl(scalable, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
int a;
struct refl_item *it = get_item(obs, refidx);
const signed int h = it->h;
const signed int k = it->k;
const signed int l = it->l;
signed int h, k, l;
get_indices(refl, &h, &k, &l);
/* For this reflection, calculate all the possible
* values of uha and vha */
......@@ -290,7 +292,15 @@ static double iterate_scale(struct image *images, int n,
if ( isnan(Ih) ) Ih = 0.0;
} else {
/* Look up by asymmetric indices */
Ih = lookup_intensity(reference, h, k, l);
Reflection *r = find_refl(reference, h, k, l);
if ( r == NULL ) {
ERROR("%3i %3i %3i isn't in the "
"reference list, so why is it "
"marked as scalable?\n", h, k, l);
Ih = 0.0;
} else {
Ih = get_intensity(r);
}
}
rha = vha - image_a->osf * uha * Ih;
......@@ -310,7 +320,9 @@ static double iterate_scale(struct image *images, int n,
if ( (reference != NULL) && (a!=b) ) continue;
/* Zero if no common reflections */
if ( cref[a+n*b] != 0 ) continue;
if ( (reference == NULL) && cref[a+n*b] != 0 ) {
continue;
}
uhb = uha_arr[b];
vhb = vha_arr[b];
......@@ -335,7 +347,6 @@ static double iterate_scale(struct image *images, int n,
gsl_vector_set(v, a, vval+vc);
}
}
free(uh_arr);
......@@ -379,27 +390,28 @@ static double iterate_scale(struct image *images, int n,
}
static RefList *lsq_intensities(struct image *images, int n,
ReflItemList *obs, const char *sym)
static void lsq_intensities(struct image *images, int n, RefList *full)
{
RefList *full;
int i;
full = reflist_new();
for ( i=0; i<num_items(obs); i++ ) {
Reflection *refl;
RefListIterator *iter;
struct refl_item *it = get_item(obs, i);
for ( refl = first_refl(full, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
double num = 0.0;
double den = 0.0;
int m;
int redundancy = 0;
Reflection *new;
signed int h, k, l;
get_indices(refl, &h, &k, &l);
/* For each frame */
for ( m=0; m<n; m++ ) {
double G;
Reflection *refl;
Reflection *f;
/* Don't scale intensities from this image if
* post refinement failed on the last step. */
......@@ -407,18 +419,16 @@ static RefList *lsq_intensities(struct image *images, int n,
G = images[m].osf;
for ( refl = find_refl(images[m].reflections,
it->h, it->k, it->l);
refl != NULL;
refl = next_found_refl(refl) )
for ( f = find_refl(images[m].reflections, h, k, l);
f != NULL; f = next_found_refl(refl) )
{
double p;
if ( !get_scalable(refl) ) continue;
p = get_partiality(refl);
if ( !get_scalable(f) ) continue;
p = get_partiality(f);
num += get_intensity(refl) * p * G;
num += get_intensity(f) * p * G;
den += pow(p, 2.0) * pow(G, 2.0);
redundancy++;
......@@ -426,19 +436,17 @@ static RefList *lsq_intensities(struct image *images, int n,
}
if ( !isnan(num/den) ) {
new = add_refl(full, it->h, it->k, it->l);
set_int(new, num/den);
set_redundancy(new, redundancy);
} else {
ERROR("Couldn't calculate LSQ full intensity for"
"%3i %3i %3i\n", it->h, it->k, it->l);
/* Doom is probably impending */
set_int(refl, num/den);
//STATUS("%3i %3i %3i %i %i\n", h, k, l,
// redundancy, get_redundancy(refl));
if ( get_redundancy(refl) != redundancy ) {
ERROR("Didn't find all the expected parts for"
" %3i %3i %3i\n", h, k, l);
}
}
return full;
}
......@@ -459,23 +467,78 @@ static void normalise_osfs(struct image *images, int n)
}
static RefList *guess_scaled_reflections(struct image *images, int n)
{
RefList *scalable;
int i;
scalable = reflist_new();
for ( i=0; i<n; i++ ) {
Reflection *refl;
RefListIterator *iter;
/* Don't scale intensities from this image if
* post refinement failed on the last step. */
if ( images[i].pr_dud ) continue;
for ( refl = first_refl(images[i].reflections, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
if ( get_scalable(refl) ) {
signed int h, k, l;
Reflection *f;
int red;
get_indices(refl, &h, &k, &l);
f = find_refl(scalable, h, k, l);
if ( f == NULL ) {
f = add_refl(scalable, h, k, l);
}
red = get_redundancy(f);
set_redundancy(f, red+1);
}
}
}
return scalable;
}
static void show_scale_factors(struct image *images, int n)
{
int i;
for ( i=0; i<n; i++ ) {
STATUS("Image %4i: scale factor %5.2f\n", i, images[i].osf);
}
}
/* Scale the stack of images */
RefList *scale_intensities(struct image *images, int n, const char *sym,
ReflItemList *obs, char *cref, double *reference)
RefList *scale_intensities(struct image *images, int n, RefList *reference)
{
int m;
RefList *full;
int i;
double max_shift;
char *cref;
/* Start with all scale factors equal */
for ( m=0; m<n; m++ ) images[m].osf = 1.0;
cref = find_common_reflections(images, n);
full = guess_scaled_reflections(images, n);
/* Iterate */
i = 0;
do {
max_shift = iterate_scale(images, n, obs, sym, cref, reference);
max_shift = iterate_scale(images, n, full, cref, reference);
STATUS("Scaling iteration %2i: max shift = %5.2f\n",
i+1, max_shift);
i++;
......@@ -483,6 +546,10 @@ RefList *scale_intensities(struct image *images, int n, const char *sym,
} while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
full = lsq_intensities(images, n, obs, sym);
free(cref);
show_scale_factors(images, n);
lsq_intensities(images, n, full);
return full;
}
src/hrs-scaling.h
View file @ 7a7cf035
......@@ -20,10 +20,8 @@
#include "image.h"
extern RefList *scale_intensities(struct image *images, int n, const char *sym,
ReflItemList *obs, char *cref,
double *reference);
extern RefList *scale_intensities(struct image *images, int n,
RefList *reference);
extern char *find_common_reflections(struct image *images, int n);
#endif /* HRS_SCALING_H */
src/partialator.c
View file @ 7a7cf035
......@@ -66,8 +66,6 @@ static void show_help(const char *s)
struct refine_args
{
const char *sym;
ReflItemList *obs;
RefList *full;
struct image *image;
FILE *graph;
......@@ -91,7 +89,7 @@ static void refine_image(void *task, int id)
struct image *image = pargs->image;
image->id = id;
pr_refine(image, pargs->full, pargs->sym);
pr_refine(image, pargs->full);
}
......@@ -123,15 +121,13 @@ static void done_image(void *vqargs, void *task)
static void refine_all(struct image *images, int n_total_patterns,
struct detector *det, const char *sym,
ReflItemList *obs, RefList *full, int nthreads,
struct detector *det,
RefList *full, int nthreads,
FILE *graph, FILE *pgraph)
{
struct refine_args task_defaults;
struct queue_args qargs;
task_defaults.sym = sym;
task_defaults.obs = obs;
task_defaults.full = full;
task_defaults.image = NULL;
task_defaults.graph = graph;
......@@ -152,7 +148,7 @@ static void refine_all(struct image *images, int n_total_patterns,
/* Decide which reflections can be scaled */
static int select_scalable_reflections(RefList *list, ReflItemList *sc_l)
static int select_scalable_reflections(RefList *list, RefList *reference)
{
Reflection *refl;
RefListIterator *iter;
......@@ -162,29 +158,24 @@ static int select_scalable_reflections(RefList *list, ReflItemList *sc_l)
refl != NULL;
refl = next_refl(refl, iter) ) {
int scalable = 1;
int sc = 1;
double v;
if ( get_partiality(refl) < 0.1 ) scalable = 0;
if ( get_partiality(refl) < 0.1 ) sc = 0;
v = fabs(get_intensity(refl));
if ( v < 0.1 ) scalable = 0;
set_scalable(refl, scalable);
if ( scalable ) {
if ( v < 0.1 ) sc = 0;
/* If we are scaling against a reference set, we additionally
* require that this reflection is in the reference list. */
if ( reference != NULL ) {
signed int h, k, l;
nobs++;
/* Add (asymmetric) indices to list */
get_indices(refl, &h, &k, &l);
if ( !find_item(sc_l, h, k, l) ) {
add_item(sc_l, h, k, l);
}
if ( find_refl(reference, h, k, l) == NULL ) sc = 0;
}
set_scalable(refl, sc);
if ( sc ) nobs++;
}
return nobs;
......@@ -192,7 +183,7 @@ static int select_scalable_reflections(RefList *list, ReflItemList *sc_l)
static void select_reflections_for_refinement(struct image *images, int n,
ReflItemList *scalable)
RefList *full, int have_reference)
{
int i;
......@@ -200,6 +191,11 @@ static void select_reflections_for_refinement(struct image *images, int n,
Reflection *refl;
RefListIterator *iter;
int n_acc = 0;
int n_nomatch = 0;
int n_noscale = 0;
int n_fewmatch = 0;
int n_ref = 0;
for ( refl = first_refl(images[i].reflections, &iter);
refl != NULL;
......@@ -208,14 +204,47 @@ static void select_reflections_for_refinement(struct image *images, int n,
signed int h, k, l;
int sc;
n_ref++;
/* We require that the reflection itself is scalable
* (i.e. sensible partiality and intensity) and that
* the "full" estimate of this reflection is made from
* at least two parts. */
get_indices(refl, &h, &k, &l);
sc = get_scalable(refl);
if ( !sc ) {
n_noscale++;
} else {
Reflection *f = find_refl(full, h, k, l);
if ( f != NULL ) {
int r = get_redundancy(f);
if ( (r >= 2) || have_reference ) {
set_refinable(refl, 1);
n_acc++;
} else {
n_fewmatch++;
}
} else {
n_nomatch++;
}
if ( sc && find_item(scalable, h, k, l) ) {
set_refinable(refl, 1);
}
}
STATUS("Image %4i: %i guide reflections accepted "
"(%i not scalable, %i few matches, %i total)\n",
i, n_acc, n_noscale, n_fewmatch, n_ref);
/* This would be a silly situation, since there must be a match
* if THIS pattern has a scalable part of the reflection! */
assert(n_nomatch == 0);
}
}
......@@ -230,7 +259,6 @@ int main(int argc, char *argv[])
FILE *fh;
int nthreads = 1;
struct detector *det;
ReflItemList *scalable;
int i;
int n_total_patterns;
struct image *images;
......@@ -240,13 +268,12 @@ int main(int argc, char *argv[])
int n_found = 0;
int n_expected = 0;
int n_notfound = 0;
char *cref;
int n_usable_patterns = 0;
int nobs;
char *reference_file = NULL;
double *reference = NULL;
RefList *reference_list = NULL;
RefList *reference = NULL;
int n_dud;
int have_reference = 0;
/* Long options */
const struct option longopts[] = {
......@@ -359,9 +386,9 @@ int main(int argc, char *argv[])
list = read_reflections(reference_file);
free(reference_file);
if ( list == NULL ) return 1;
reference_list = asymmetric_indices(list, sym);
reference = asymmetric_indices(list, sym);
reflist_free(list);
reference = intensities_from_list(reference_list);
have_reference = 1;
}
......@@ -382,7 +409,6 @@ int main(int argc, char *argv[])
/* Fill in what we know about the images so far */
rewind(fh);
scalable = new_items();
nobs = 0;
for ( i=0; i<n_total_patterns; i++ ) {
......@@ -427,9 +453,10 @@ int main(int argc, char *argv[])
cur->reflections = as;
predict_corresponding_reflections(cur, sym, &n_expected,
&n_found, &n_notfound);
&n_found, &n_notfound);
nobs += select_scalable_reflections(cur->reflections, scalable);
nobs += select_scalable_reflections(cur->reflections,
reference);
progress_bar(i, n_total_patterns-1, "Loading pattern data");
n_usable_patterns++;
......@@ -438,51 +465,13 @@ int main(int argc, char *argv[])
fclose(fh);
STATUS("Found %5.2f%% of the expected peaks (missed %i of %i).\n",
100.0 * (double)n_found / n_expected, n_notfound, n_expected);
STATUS("Mean measurements per scalable unique reflection: %5.2f\n",
(double)nobs / num_items(scalable));
cref = find_common_reflections(images, n_usable_patterns);
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
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);
f = find_refl(full, it->h, it->k, it->l);
if ( f == NULL ) {
ERROR("%3i %3i %3i was designated scalable, but no"
" full intensity was recorded.\n",
it->h, it->k, it->l);
}
}
full = scale_intensities(images, n_usable_patterns, reference);
for ( i=0; i<n_usable_patterns; 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;
if ( !get_scalable(refl) ) continue;
get_indices(refl, &h, &k, &l);
if ( find_item(scalable, h, k, l) == 0 ) {
ERROR("%3i %3i %3i in image %i is scalable"
" but is not in the list of scalable"
" reflections.\n", h, k, l, i);
}
}
}
select_reflections_for_refinement(images, n_usable_patterns, full,