diff options
-rw-r--r-- | src/likelihood.c | 9 | ||||
-rw-r--r-- | src/likelihood.h | 7 | ||||
-rw-r--r-- | src/process_hkl.c | 37 | ||||
-rw-r--r-- | src/utils.c | 90 | ||||
-rw-r--r-- | src/utils.h | 22 |
5 files changed, 144 insertions, 21 deletions
diff --git a/src/likelihood.c b/src/likelihood.c index 95f005c7..98ffcba0 100644 --- a/src/likelihood.c +++ b/src/likelihood.c @@ -19,17 +19,20 @@ void detwin_intensities(const double *model, double *new_pattern, const unsigned int *model_counts, - unsigned int *new_counts) + ReflItemList *items) { /* Placeholder... */ } void scale_intensities(const double *model, double *new_pattern, const unsigned int *model_counts, - unsigned int *new_counts, double f0, int f0_valid) + ReflItemList *items, double f0, int f0_valid) { double s; unsigned int i; + unsigned int *new_counts; + + new_counts = items_to_counts(items); s = stat_scale_intensity(model, model_counts, new_pattern, new_counts); if ( f0_valid ) printf("%f %f\n", s, f0); @@ -41,4 +44,6 @@ void scale_intensities(const double *model, double *new_pattern, for ( i=0; i<LIST_SIZE; i++ ) { new_counts[i] *= s; } + + free(new_counts); } diff --git a/src/likelihood.h b/src/likelihood.h index 0a2bb407..ad6585eb 100644 --- a/src/likelihood.h +++ b/src/likelihood.h @@ -17,13 +17,16 @@ #define LIKELIHOOD_H +#include "utils.h" + + extern void detwin_intensities(const double *model, double *new_pattern, const unsigned int *model_counts, - unsigned int *new_counts); + ReflItemList *items); extern void scale_intensities(const double *model, double *new_pattern, const unsigned int *model_counts, - unsigned int *new_counts, double f0, + ReflItemList *items, double f0, int f0_valid); diff --git a/src/process_hkl.c b/src/process_hkl.c index 77d5056a..deeb77f5 100644 --- a/src/process_hkl.c +++ b/src/process_hkl.c @@ -188,17 +188,21 @@ static void process_reflections(double *ref, unsigned int *counts, static void merge_pattern(double *model, const double *new, unsigned int *model_counts, - const unsigned int *counts, int mo, int sum) + ReflItemList *items, int mo, int sum) { - signed int h, k, l; + int i; - for ( l=-INDMAX; l<INDMAX; l++ ) { - for ( k=-INDMAX; k<INDMAX; k++ ) { - for ( h=-INDMAX; h<INDMAX; h++ ) { + for ( i=0; i<num_items(items); i++ ) { double intensity; + signed int h, k, l; + struct refl_item *item; - if ( lookup_count(counts, h, k, l) == 0 ) continue; + item = get_item(items, i); + + h = item->h; + k = item->k; + l = item->l; intensity = lookup_intensity(new, h, k, l); @@ -217,8 +221,6 @@ static void merge_pattern(double *model, const double *new, } } - } - } } @@ -243,7 +245,6 @@ int main(int argc, char *argv[]) int config_scale = 0; char *intfile = NULL; double *new_pattern = NULL; - unsigned int *new_counts = NULL; unsigned int n_total_patterns; unsigned int *truecounts = NULL; char *sym = NULL; @@ -251,6 +252,7 @@ int main(int argc, char *argv[]) float f0; int f0_valid; int n_nof0 = 0; + ReflItemList *items; /* Long options */ const struct option longopts[] = { @@ -344,7 +346,7 @@ int main(int argc, char *argv[]) model_counts = new_list_count(); cell = load_cell_from_pdb(pdb); new_pattern = new_list_intensity(); - new_counts = new_list_count(); + items = new_items(); if ( strcmp(filename, "-") == 0 ) { fh = stdin; @@ -393,7 +395,7 @@ int main(int argc, char *argv[]) /* Detwin before scaling */ if ( config_detwin ) { detwin_intensities(model, new_pattern, - model_counts, new_counts); + model_counts, items); } /* Assume a default I0 if we don't have one by now */ @@ -405,13 +407,13 @@ int main(int argc, char *argv[]) /* Scale if requested */ if ( config_scale ) { scale_intensities(model, new_pattern, - model_counts, new_counts, f0, + model_counts, items, f0, f0_valid); } /* Start of second or later pattern */ merge_pattern(model, new_pattern, model_counts, - new_counts, config_maxonly, config_sum); + items, config_maxonly, config_sum); if ( (trueref != NULL) && config_every && (n_patterns % config_every == 0) ) { @@ -424,9 +426,8 @@ int main(int argc, char *argv[]) if ( n_patterns == config_stopafter ) break; - zero_list_count(new_counts); - n_patterns++; + clear_items(items); progress_bar(n_patterns, n_total_patterns, "Merging"); @@ -449,12 +450,12 @@ int main(int argc, char *argv[]) if ( (h==0) && (k==0) && (l==0) ) continue; - if ( lookup_count(new_counts, h, k, l) != 0 ) { + if ( find_item(items, h, k, l) != 0 ) { ERROR("More than one measurement for %i %i %i in" " pattern number %i\n", h, k, l, n_patterns); } set_intensity(new_pattern, h, k, l, intensity); - set_count(new_counts, h, k, l, 1); + add_item(items, h, k, l); } while ( rval != NULL ); @@ -480,5 +481,7 @@ int main(int argc, char *argv[]) STATUS("There were %u patterns.\n", n_patterns); STATUS("%i had no f0 valid value.\n", n_nof0); + delete_items(items); + return 0; } diff --git a/src/utils.c b/src/utils.c index e3b8b667..de820cf3 100644 --- a/src/utils.c +++ b/src/utils.c @@ -283,3 +283,93 @@ int assplode(const char *a, const char *delims, char ***pbits, *pbits = bits; return n; } + + +struct _reflitemlist { + struct refl_item *items; + int n_items; + int max_items; +}; + +void clear_items(ReflItemList *items) +{ + items->n_items = 0; +} + +static void alloc_items(ReflItemList *items) +{ + items->items = realloc(items->items, + items->max_items*sizeof(struct refl_item)); +} + +ReflItemList *new_items() +{ + ReflItemList *new; + new = malloc(sizeof(ReflItemList)); + new->max_items = 1024; + new->n_items = 0; + new->items = NULL; + alloc_items(new); + return new; +} + +void delete_items(ReflItemList *items) +{ + if ( items->items != NULL ) free(items->items); + free(items); +} + +void add_item(ReflItemList *items, + signed int h, signed int k, signed int l) +{ + if ( items->n_items == items->max_items ) { + items->max_items += 1024; + alloc_items(items); + } + + items->items[items->n_items].h = h; + items->items[items->n_items].k = k; + items->items[items->n_items].l = l; + items->n_items++; +} + +int find_item(ReflItemList *items, + signed int h, signed int k, signed int l) +{ + int i; + + for ( i=0; i<items->n_items; i++ ) { + if ( items->items[i].h != h ) continue; + if ( items->items[i].k != k ) continue; + if ( items->items[i].l != l ) continue; + return 1; + } + return 0; +} + +struct refl_item *get_item(ReflItemList *items, int i) +{ + if ( i >= items->n_items ) return NULL; + return &items->items[i]; +} + +int num_items(ReflItemList *items) +{ + return items->n_items; +} + +unsigned int *items_to_counts(ReflItemList *items) +{ + unsigned int *c; + int i; + + c = new_list_count(); + + for ( i=0; i<num_items(items); i++ ) { + struct refl_item *r; + r = get_item(items, i); + set_count(c, r->h, r->k, r->l, 1); + } + + return c; +} diff --git a/src/utils.h b/src/utils.h index b01725dd..14ef45a7 100644 --- a/src/utils.h +++ b/src/utils.h @@ -176,6 +176,28 @@ static inline double angle_between(double x1, double y1, double z1, #include "list_tmp.h" +/* ----------- Reflection lists indexed by sequence (not indices) ----------- */ + +typedef struct _reflitemlist ReflItemList; /* Opaque */ + +struct refl_item { + signed int h; + signed int k; + signed int l; +}; + +extern void clear_items(ReflItemList *items); +extern ReflItemList *new_items(void); +extern void delete_items(ReflItemList *items); +extern void add_item(ReflItemList *items, + signed int h, signed int k, signed int l); +extern int find_item(ReflItemList *items, + signed int h, signed int k, signed int l); +extern struct refl_item *get_item(ReflItemList *items, int i); +extern int num_items(ReflItemList *items); +extern unsigned int *items_to_counts(ReflItemList *items); + + /* ------------------------------ Message macros ---------------------------- */ #define ERROR(...) fprintf(stderr, __VA_ARGS__) |