From e980ed54dc29e025587aba47390727c500aec8f1 Mon Sep 17 00:00:00 2001 From: Thomas White Date: Tue, 8 Feb 2011 19:10:27 +0100 Subject: Work on making iteration work --- src/geometry.c | 5 +- src/hrs-scaling.c | 3 + src/partialator.c | 12 ++-- src/peaks.c | 5 +- src/post-refinement.c | 5 +- src/reflist.c | 170 ++++++++++++++++++++++++++++++++++++++++---------- src/reflist.h | 5 +- src/templates.c | 10 +-- tests/list_check.c | 49 ++++++++++++++- 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; hh, 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; mreflections == 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