aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/geometry.c5
-rw-r--r--src/hrs-scaling.c3
-rw-r--r--src/partialator.c12
-rw-r--r--src/peaks.c5
-rw-r--r--src/post-refinement.c5
-rw-r--r--src/reflist.c170
-rw-r--r--src/reflist.h5
-rw-r--r--src/templates.c10
-rw-r--r--tests/list_check.c49
9 files changed, 212 insertions, 52 deletions
diff --git a/src/geometry.c b/src/geometry.c
index 9f9b953c..1cf839c2 100644
--- a/src/geometry.c
+++ b/src/geometry.c
@@ -263,10 +263,11 @@ double integrate_all(struct image *image, RefList *reflections)
{
double itot = 0.0;
Reflection *refl;
+ RefListIterator *iter;
- for ( refl = first_refl(reflections);
+ for ( refl = first_refl(reflections, &iter);
refl != NULL;
- refl = next_refl(refl) ) {
+ refl = next_refl(refl, iter) ) {
float x, y, intensity;
double xp, yp;
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index d5d96ffc..9f1ada75 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -52,6 +52,7 @@ static void s_uhavha(signed int hat, signed int kat, signed int lat,
if ( !get_scalable(refl) ) continue;
ic = get_intensity(refl) / get_partiality(refl);
+ STATUS("%f / %f = %f\n", get_intensity(refl), get_partiality(refl), ic);
sigi = sqrt(fabs(ic));
uha_val += 1.0 / pow(sigi, 2.0);
@@ -117,6 +118,7 @@ static double iterate_scale(struct image *images, int n,
uh_arr = new_list_intensity();
vh_arr = new_list_intensity();
+ STATUS("%i\n", n_ref);
for ( h=0; h<n_ref; h++ ) {
double uh, vh;
@@ -124,6 +126,7 @@ static double iterate_scale(struct image *images, int n,
uh = s_uh(images, n, it->h, it->k, it->l, sym);
vh = s_vh(images, n, it->h, it->k, it->l, sym);
+ STATUS("%f %f\n", uh, vh);
set_intensity(uh_arr, it->h, it->k, it->l, uh);
set_intensity(vh_arr, it->h, it->k, it->l, vh);
diff --git a/src/partialator.c b/src/partialator.c
index 37c8fbd6..3251f7e8 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -116,10 +116,11 @@ static void select_scalable_reflections(struct image *images, int n)
for ( m=0; m<n; m++ ) {
Reflection *refl;
+ RefListIterator *iter;
- for ( refl = first_refl(images[m].reflections);
+ for ( refl = first_refl(images[m].reflections, &iter);
refl != NULL;
- refl = next_refl(refl) ) {
+ refl = next_refl(refl, iter) ) {
int scalable = 1;
double v;
@@ -275,6 +276,7 @@ int main(int argc, char *argv[])
RefList *peaks;
RefList *transfer;
Reflection *refl;
+ RefListIterator *iter;
if ( find_chunk(fh, &cell, &filename) == 1 ) {
ERROR("Couldn't get all of the filenames and cells"
@@ -327,9 +329,9 @@ int main(int argc, char *argv[])
transfer = find_intersections(&images[i], cell, 0);
images[i].reflections = reflist_new();
- for ( refl = first_refl(transfer);
+ for ( refl = first_refl(transfer, &iter);
refl != NULL;
- refl = next_refl(refl) ) {
+ refl = next_refl(refl, iter) ) {
Reflection *peak;
Reflection *new;
@@ -338,6 +340,8 @@ int main(int argc, char *argv[])
int clamp1, clamp2;
get_indices(refl, &h, &k, &l);
+ STATUS("%3i %3i %3i\n", h, k, l);
+
peak = find_refl(peaks, h, k, l);
if ( peak == NULL ) {
if ( (h==0) && (k==0) && (l==0) ) continue;
diff --git a/src/peaks.c b/src/peaks.c
index 06c5ee35..4ae0a70b 100644
--- a/src/peaks.c
+++ b/src/peaks.c
@@ -669,6 +669,7 @@ void output_intensities(struct image *image, UnitCell *cell,
double bsx, bsy, bsz;
double csx, csy, csz;
Reflection *refl;
+ RefListIterator *iter;
if ( image->reflections == NULL ) {
find_projected_peaks(image, cell, circular_domain, domain_r);
@@ -683,9 +684,9 @@ void output_intensities(struct image *image, UnitCell *cell,
&bsx, &bsy, &bsz,
&csx, &csy, &csz);
- for ( refl = first_refl(image->reflections);
+ for ( refl = first_refl(image->reflections, &iter);
refl != NULL;
- refl = next_refl(refl) ) {
+ refl = next_refl(refl, iter) ) {
float x, y, intensity;
double d;
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 4fa7ad47..0bd26f76 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -257,6 +257,7 @@ static double pr_iterate(struct image *image, const double *i_full,
gsl_vector *shifts;
int param;
Reflection *refl;
+ RefListIterator *iter;
RefList *reflections;
double max_shift;
@@ -266,9 +267,9 @@ static double pr_iterate(struct image *image, const double *i_full,
v = gsl_vector_calloc(NUM_PARAMS);
/* Construct the equations, one per reflection in this image */
- for ( refl = first_refl(reflections);
+ for ( refl = first_refl(reflections, &iter);
refl != NULL;
- refl = next_refl(refl) ) {
+ refl = next_refl(refl, iter) ) {
signed int hind, kind, lind;
signed int ha, ka, la;
diff --git a/src/reflist.c b/src/reflist.c
index e543c20d..4b08b3af 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -24,6 +24,7 @@ struct _reflection {
struct _reflection *parent; /* Parent node */
struct _reflection *next; /* Another reflection with the same
* indices, or NULL */
+ struct _reflection *prev;
signed int h;
signed int k;
@@ -71,6 +72,7 @@ static Reflection *new_node(unsigned int serial)
new = calloc(1, sizeof(struct _reflection));
new->serial = serial;
new->next = NULL;
+ new->prev = NULL;
new->child[0] = NULL;
new->child[1] = NULL;
@@ -110,24 +112,36 @@ void reflist_free(RefList *list)
/********************************** Search ************************************/
-static Reflection *recursive_search(Reflection *refl, unsigned int search)
+/* Return the first reflection in 'list' with the given indices, or NULL */
+Reflection *find_refl(RefList *list, INDICES)
{
- int i;
+ unsigned int search = SERIAL(h, k, l);
+ Reflection *refl = list->head->child[0];
- /* Hit the bottom of the tree? */
- if ( refl == NULL ) return NULL;
+ while ( refl != NULL ) {
- /* Is this the correct reflection? */
- if ( refl->serial == search ) return refl;
+ if ( search < refl->serial ) {
- for ( i=0; i<2; i++ ) {
+ if ( refl->child[0] != NULL ) {
+ refl = refl->child[0];
+ } else {
+ /* Hit the bottom of the tree */
+ return NULL;
+ }
- if ( search < refl->serial ) {
- return recursive_search(refl->child[0], search);
- }
+ } else if ( search > refl->serial ) {
+
+ if ( refl->child[1] != NULL ) {
+ refl = refl->child[1];
+ } else {
+ /* Hit the bottom of the tree */
+ return NULL;
+ }
+
+ } else {
+
+ return refl;
- if ( search >= refl->serial ) {
- return recursive_search(refl->child[1], search);
}
}
@@ -136,14 +150,6 @@ static Reflection *recursive_search(Reflection *refl, unsigned int search)
}
-/* Return the first reflection in 'list' with the given indices, or NULL */
-Reflection *find_refl(RefList *list, INDICES)
-{
- unsigned int search = SERIAL(h, k, l);
- return recursive_search(list->head, search);
-}
-
-
/* Find the next reflection in 'refl's list with the same indices, or NULL */
Reflection *next_found_refl(Reflection *refl)
{
@@ -247,6 +253,7 @@ static void insert_node(Reflection *head, Reflection *new)
while ( refl != NULL ) {
if ( new->serial < refl->serial ) {
+
if ( refl->child[0] != NULL ) {
refl = refl->child[0];
} else {
@@ -254,7 +261,9 @@ static void insert_node(Reflection *head, Reflection *new)
new->parent = refl;
return;
}
- } else {
+
+ } else if ( new->serial > refl->serial ) {
+
if ( refl->child[1] != NULL ) {
refl = refl->child[1];
} else {
@@ -262,8 +271,20 @@ static void insert_node(Reflection *head, Reflection *new)
new->parent = refl;
return;
}
- }
+ } else {
+
+ /* New reflection is identical to a previous one */
+ do {
+ if ( refl->next == NULL ) {
+ refl->next = new;
+ new->prev = refl;
+ return;
+ }
+ refl = refl->next;
+ } while ( 1 );
+
+ }
}
}
@@ -294,6 +315,33 @@ void delete_refl(Reflection *refl)
int i;
Reflection **parent_pos = NULL;
+ /* Is this a duplicate, and not the first one? */
+ if ( refl->prev != NULL ) {
+ refl->prev->next = refl->next;
+ free(refl);
+ return;
+ }
+
+ /* Is this the first duplicate of many? */
+ if ( refl->next != NULL ) {
+
+ refl->next->parent = refl->parent;
+ refl->next->child[0] = refl->child[0];
+ refl->next->child[1] = refl->child[1];
+ refl->next->prev = NULL;
+
+ for ( i=0; i<2; i++ ) {
+ if ( refl->parent->child[i] == refl ) {
+ refl->parent->child[i] = refl->next;
+ }
+ }
+
+ free(refl);
+
+ return;
+
+ }
+
/* Remove parent's reference */
for ( i=0; i<2; i++ ) {
if ( refl->parent->child[i] == refl ) {
@@ -344,24 +392,80 @@ void delete_refl(Reflection *refl)
/********************************* Iteration **********************************/
-Reflection *first_refl(RefList *list)
+
+struct _reflistiterator {
+
+ int stack_size;
+ int stack_ptr;
+ Reflection **stack;
+
+};
+
+Reflection *first_refl(RefList *list, RefListIterator **piter)
{
- return list->head->child[0];
+ RefListIterator *iter;
+
+ iter = malloc(sizeof(struct _reflistiterator));
+ iter->stack_size = 32;
+ iter->stack = malloc(iter->stack_size*sizeof(Reflection *));
+ iter->stack_ptr = 0;
+ *piter = iter;
+
+ Reflection *refl = list->head->child[0];
+
+ do {
+
+ if ( refl != NULL ) {
+ iter->stack[iter->stack_ptr++] = refl;
+ if ( iter->stack_ptr == iter->stack_size ) {
+ iter->stack_size += 32;
+ iter->stack = realloc(iter->stack,
+ iter->stack_size*sizeof(Reflection *));
+ }
+ refl = refl->child[0];
+ continue;
+ }
+
+ if ( iter->stack_ptr == 0 ) {
+ free(iter->stack);
+ free(iter);
+ return NULL;
+ }
+
+ refl = iter->stack[iter->stack_ptr--];
+
+ return refl;
+
+ } while ( 1 );
}
-Reflection *next_refl(Reflection *refl)
+Reflection *next_refl(Reflection *refl, RefListIterator *iter)
{
- /* Does a left child exist? */
- if ( refl->child[0] != NULL ) return refl->child[0];
+ do {
- /* Otherwise move up the tree to find the next right child */
- while ( refl->child[1] != NULL ) {
- refl = refl->parent;
- if ( refl == NULL ) return NULL;
- }
+ refl = refl->child[1];;
+
+ if ( refl != NULL ) {
+ iter->stack[iter->stack_ptr++] = refl;
+ if ( iter->stack_ptr == iter->stack_size ) {
+ iter->stack_size += 32;
+ iter->stack = realloc(iter->stack,
+ iter->stack_size*sizeof(Reflection *));
+ }
+ refl = refl->child[0];
+ continue;
+ }
+ if ( iter->stack_ptr == 0 ) {
+ free(iter->stack);
+ free(iter);
+ return NULL;
+ }
+ refl = iter->stack[iter->stack_ptr--];
+
+ return refl;
- return refl->child[1];
+ } while ( 1 );
}
diff --git a/src/reflist.h b/src/reflist.h
index 78aadf2a..e25a16ad 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -19,6 +19,7 @@
typedef struct _reflist RefList;
typedef struct _reflection Reflection;
+typedef struct _reflistiterator RefListIterator;
#define INDICES signed int h, signed int k, signed int l
@@ -56,8 +57,8 @@ extern Reflection *add_refl(RefList *list, INDICES);
extern void delete_refl(Reflection *refl);
/* Iteration */
-extern Reflection *first_refl(RefList *list);
-extern Reflection *next_refl(Reflection *refl);
+extern Reflection *first_refl(RefList *list, RefListIterator **iterator);
+extern Reflection *next_refl(Reflection *refl, RefListIterator *iter);
/* Voodoo */
extern void optimise_reflist(RefList *list);
diff --git a/src/templates.c b/src/templates.c
index 9abbcb17..79e60aec 100644
--- a/src/templates.c
+++ b/src/templates.c
@@ -232,13 +232,14 @@ static double integrate_all_rot(struct image *image, RefList *reflections,
double cosr, sinr;
int num_int = 0;
Reflection *refl;
+ RefListIterator *iter;
cosr = cos(rot);
sinr = sin(rot);
- for ( refl = first_refl(reflections);
+ for ( refl = first_refl(reflections, &iter);
refl != NULL;
- refl = next_refl(refl) ) {
+ refl = next_refl(refl, iter) ) {
float xp, yp;
double x, y;
@@ -266,14 +267,15 @@ static double mean_distance(struct image *image, RefList *reflections,
double cosr, sinr;
int num_dist = 0;
Reflection *refl;
+ RefListIterator *iter;
cosr = cos(rot);
sinr = sin(rot);
/* For each template peak */
- for ( refl = first_refl(reflections);
+ for ( refl = first_refl(reflections, &iter);
refl != NULL;
- refl = next_refl(refl) ) {
+ refl = next_refl(refl, iter) ) {
float xp, yp;
int j;
diff --git a/tests/list_check.c b/tests/list_check.c
index 00e27e00..ffae8d12 100644
--- a/tests/list_check.c
+++ b/tests/list_check.c
@@ -27,6 +27,7 @@ struct refltemp {
signed int l;
int del;
int dup;
+ int found;
};
#define RANDOM_INDEX (128*random()/RAND_MAX - 256*random()/RAND_MAX)
@@ -38,6 +39,8 @@ static int test_lists(int num_items)
RefList *list;
int i;
signed int h, k, l;
+ Reflection *refl;
+ RefListIterator *iter;
check = malloc(num_items * sizeof(struct refltemp));
list = reflist_new();
@@ -53,11 +56,11 @@ static int test_lists(int num_items)
int j;
int duplicate = 0;
- if ( random() > RAND_MAX/2 ) {
+ //if ( random() > RAND_MAX/2 ) {
h = RANDOM_INDEX;
k = RANDOM_INDEX;
l = RANDOM_INDEX;
- } /* else use the same as last time */
+ //} /* else use the same as last time */
/* Count the number of times this reflection appeared before */
for ( j=0; j<i; j++ ) {
@@ -83,9 +86,49 @@ static int test_lists(int num_items)
check[i].l = l;
check[i].del = 0;
check[i].dup = duplicate;
+ check[i].found = 0;
}
+ /* Iterate over the list and check we find everything */
+ for ( refl = first_refl(list, &iter);
+ refl != NULL;
+ refl = next_refl(refl, iter) ) {
+
+ signed int h, k, l;
+
+ get_indices(refl, &h, &k, &l);
+ printf("%3i %3i %3i\n", h, k, l);
+
+ for ( i=0; i<num_items; i++ ) {
+ if ( (check[i].h == h)
+ && (check[i].k == k)
+ && (check[i].l == l)
+ && (check[i].found == 0) ) {
+ check[i].found = 1;
+ }
+ }
+
+ }
+ for ( i=0; i<num_items; i++ ) {
+ if ( check[i].found == 0 ) {
+
+ Reflection *test;
+
+ fprintf(stderr, "Iteration didn't find %3i %3i %3i %i\n",
+ check[i].h, check[i].k, check[i].l, i);
+ test = find_refl(list, check[i].h, check[i].k,
+ check[i].l);
+ if ( test == NULL ) {
+ fprintf(stderr, "Not in list\n");
+ } else {
+ fprintf(stderr, "But found in list.\n");
+ }
+ return 1;
+
+ }
+ }
+
/* Check that all the reflections can be found,
* and delete the first few. */
for ( i=0; i<num_items; i++ ) {
@@ -187,7 +230,7 @@ int main(int argc, char *argv[])
printf("Running list test...");
fflush(stdout);
- for ( i=0; i<100; i++ ) {
+ for ( i=0; i<1; i++ ) {
if ( test_lists(4096*random()/RAND_MAX) ) return 1;
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 21:33:15 +0000