aboutsummaryrefslogtreecommitdiff
path: root/src
diff options
context:
space:
mode:
authorThomas White <taw@physics.org>2011-03-22 14:59:05 +0100
committerThomas White <taw@physics.org>2012-02-22 15:27:20 +0100
commit656abb22b20cf9a4a47f6062e80ffa9553d38c0f (patch)
treec5079ff1f9af0bed2ed86ff0b0594daa9e63ceff /src
parent407a74d5f4ad651170707ae762a63c6cbec47205 (diff)
partialator: Save reflection list
Diffstat (limited to 'src')
-rw-r--r--src/hrs-scaling.c24
-rw-r--r--src/hrs-scaling.h4
-rw-r--r--src/partialator.c27
-rw-r--r--src/post-refinement.c13
-rw-r--r--src/post-refinement.h2
-rw-r--r--src/reflist.c27
-rw-r--r--src/reflist.h26
7 files changed, 66 insertions, 57 deletions
diff --git a/src/hrs-scaling.c b/src/hrs-scaling.c
index e2adaaf6..03458670 100644
--- a/src/hrs-scaling.c
+++ b/src/hrs-scaling.c
@@ -311,19 +311,20 @@ static double iterate_scale(struct image *images, int n,
}
-static double *lsq_intensities(struct image *images, int n,
- ReflItemList *obs, const char *sym)
+static RefList *lsq_intensities(struct image *images, int n,
+ ReflItemList *obs, const char *sym)
{
- double *I_full;
+ RefList *full;
int i;
- I_full = new_list_intensity();
+ full = reflist_new();
for ( i=0; i<num_items(obs); i++ ) {
struct refl_item *it = get_item(obs, i);
double num = 0.0;
double den = 0.0;
int m;
+ Reflection *new;
/* For each frame */
for ( m=0; m<n; m++ ) {
@@ -351,11 +352,12 @@ static double *lsq_intensities(struct image *images, int n,
}
- set_intensity(I_full, it->h, it->k, it->l, num/den);
+ new = add_refl(full, it->h, it->k, it->l);
+ set_int(new, num/den);
}
- return I_full;
+ return full;
}
@@ -377,11 +379,11 @@ static void normalise_osfs(struct image *images, int n)
/* Scale the stack of images */
-double *scale_intensities(struct image *images, int n, const char *sym,
- ReflItemList *obs, char *cref)
+RefList *scale_intensities(struct image *images, int n, const char *sym,
+ ReflItemList *obs, char *cref)
{
int m;
- double *I_full;
+ RefList *full;
int i;
double max_shift;
@@ -399,6 +401,6 @@ double *scale_intensities(struct image *images, int n, const char *sym,
} while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
- I_full = lsq_intensities(images, n, obs, sym);
- return I_full;
+ full = lsq_intensities(images, n, obs, sym);
+ return full;
}
diff --git a/src/hrs-scaling.h b/src/hrs-scaling.h
index 04a321e3..d65129b8 100644
--- a/src/hrs-scaling.h
+++ b/src/hrs-scaling.h
@@ -20,8 +20,8 @@
#include "image.h"
-extern double *scale_intensities(struct image *images, int n, const char *sym,
- ReflItemList *obs, char *cref);
+extern RefList *scale_intensities(struct image *images, int n, const char *sym,
+ ReflItemList *obs, char *cref);
extern char *find_common_reflections(struct image *images, int n);
diff --git a/src/partialator.c b/src/partialator.c
index 553a216a..700159ee 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -64,7 +64,7 @@ struct refine_args
{
const char *sym;
ReflItemList *obs;
- double *i_full;
+ RefList *full;
struct image *image;
FILE *graph;
FILE *pgraph;
@@ -77,13 +77,13 @@ static void refine_image(int mytask, void *tasks)
struct refine_args *pargs = &all_args[mytask];
struct image *image = pargs->image;
- pr_refine(image, pargs->i_full, pargs->sym);
+ pr_refine(image, pargs->full, pargs->sym);
}
static void refine_all(struct image *images, int n_total_patterns,
struct detector *det, const char *sym,
- ReflItemList *obs, double *i_full, int nthreads,
+ ReflItemList *obs, RefList *full, int nthreads,
FILE *graph, FILE *pgraph)
{
struct refine_args *tasks;
@@ -94,7 +94,7 @@ static void refine_all(struct image *images, int n_total_patterns,
tasks[i].sym = sym;
tasks[i].obs = obs;
- tasks[i].i_full = i_full;
+ tasks[i].full = full;
tasks[i].image = &images[i];
tasks[i].graph = graph;
tasks[i].pgraph = pgraph;
@@ -157,7 +157,7 @@ int main(int argc, char *argv[])
struct image *images;
int n_iter = 10;
struct beam_params *beam = NULL;
- double *I_full;
+ RefList *full;
int n_found = 0;
int n_expected = 0;
int n_notfound = 0;
@@ -286,6 +286,9 @@ int main(int argc, char *argv[])
return 1;
}
+ /* FIXME: This leaves gaps in the array */
+ if ( images[i].indexed_cell == NULL ) continue;
+
/* Won't be needing this, if it exists */
image_feature_list_free(images[i].features);
images[i].features = NULL;
@@ -372,7 +375,7 @@ int main(int argc, char *argv[])
/* Make initial estimates */
STATUS("Performing initial scaling.\n");
select_scalable_reflections(images, n_total_patterns);
- I_full = scale_intensities(images, n_total_patterns, sym, obs, cref);
+ full = scale_intensities(images, n_total_patterns, sym, obs, cref);
/* Iterate */
for ( i=0; i<n_iter; i++ ) {
@@ -398,14 +401,14 @@ int main(int argc, char *argv[])
}
/* Refine the geometry of all patterns to get the best fit */
- refine_all(images, n_total_patterns, det, sym, obs, I_full,
+ refine_all(images, n_total_patterns, det, sym, obs, full,
nthreads, fhg, fhp);
/* Re-estimate all the full intensities */
- free(I_full);
+ reflist_free(full);
select_scalable_reflections(images, n_total_patterns);
- I_full = scale_intensities(images, n_total_patterns,
- sym, obs, cref);
+ full = scale_intensities(images, n_total_patterns,
+ sym, obs, cref);
fclose(fhg);
fclose(fhp);
@@ -418,13 +421,13 @@ int main(int argc, char *argv[])
}
/* Output results */
- //write_reflist(outfile, obs, I_full, NULL, NULL, cts, NULL);
+ write_reflist(outfile, full, images[0].indexed_cell);
/* Clean up */
for ( i=0; i<n_total_patterns; i++ ) {
reflist_free(images[i].reflections);
}
- free(I_full);
+ reflist_free(full);
delete_items(obs);
free(sym);
free(outfile);
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 0bd26f76..d5f9d3f8 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -248,8 +248,8 @@ static void apply_shift(struct image *image, int k, double shift)
}
-/* Perform one cycle of post refinement on 'image' against 'i_full' */
-static double pr_iterate(struct image *image, const double *i_full,
+/* Perform one cycle of post refinement on 'image' against 'full' */
+static double pr_iterate(struct image *image, const RefList *full,
const char *sym)
{
gsl_matrix *M;
@@ -277,6 +277,7 @@ static double pr_iterate(struct image *image, const double *i_full,
double I_partial;
int k;
double p;
+ Reflection *match;
get_indices(refl, &hind, &kind, &lind);
@@ -286,7 +287,9 @@ static double pr_iterate(struct image *image, const double *i_full,
I_partial = get_intensity(refl);
get_asymm(hind, kind, lind, &ha, &ka, &la, sym);
- I_full = lookup_intensity(i_full, ha, ka, la);
+ match = find_refl(full, ha, ka, la);
+ assert(match != NULL);
+ I_full = get_intensity(match);
p = get_partiality(refl);
delta_I = I_partial - (p * I_full / image->osf);
@@ -338,14 +341,14 @@ static double pr_iterate(struct image *image, const double *i_full,
}
-void pr_refine(struct image *image, const double *i_full, const char *sym)
+void pr_refine(struct image *image, const RefList *full, const char *sym)
{
double max_shift;
int i;
i = 0;
do {
- max_shift = pr_iterate(image, i_full, sym);
+ max_shift = pr_iterate(image, full, sym);
STATUS("Iteration %2i: max shift = %5.2f\n", i, max_shift);
i++;
} while ( (max_shift > 0.01) && (i < MAX_CYCLES) );
diff --git a/src/post-refinement.h b/src/post-refinement.h
index a10fafe4..3c8d8587 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -24,7 +24,7 @@
#include "utils.h"
-extern void pr_refine(struct image *image, const double *i_full,
+extern void pr_refine(struct image *image, const RefList *full,
const char *sym);
diff --git a/src/reflist.c b/src/reflist.c
index a81c738e..30d9f696 100644
--- a/src/reflist.c
+++ b/src/reflist.c
@@ -133,7 +133,7 @@ void reflist_free(RefList *list)
/********************************** Search ************************************/
/* Return the first reflection in 'list' with the given indices, or NULL */
-Reflection *find_refl(RefList *list, INDICES)
+Reflection *find_refl(const RefList *list, INDICES)
{
unsigned int search = SERIAL(h, k, l);
Reflection *refl = list->head->child[0];
@@ -185,20 +185,21 @@ Reflection *next_found_refl(Reflection *refl)
/********************************** Getters ***********************************/
-double get_excitation_error(Reflection *refl)
+double get_excitation_error(const Reflection *refl)
{
return refl->data.excitation_error;
}
-void get_detector_pos(Reflection *refl, double *x, double *y)
+void get_detector_pos(const Reflection *refl, double *x, double *y)
{
*x = refl->data.x;
*y = refl->data.y;
}
-void get_indices(Reflection *refl, signed int *h, signed int *k, signed int *l)
+void get_indices(const Reflection *refl,
+ signed int *h, signed int *k, signed int *l)
{
*h = refl->data.h;
*k = refl->data.k;
@@ -206,19 +207,19 @@ void get_indices(Reflection *refl, signed int *h, signed int *k, signed int *l)
}
-double get_partiality(Reflection *refl)
+double get_partiality(const Reflection *refl)
{
return refl->data.p;
}
-double get_intensity(Reflection *refl)
+double get_intensity(const Reflection *refl)
{
return refl->data.intensity;
}
-void get_partial(Reflection *refl, double *r1, double *r2, double *p,
+void get_partial(const Reflection *refl, double *r1, double *r2, double *p,
int *clamp_low, int *clamp_high)
{
*r1 = refl->data.r1;
@@ -229,31 +230,31 @@ void get_partial(Reflection *refl, double *r1, double *r2, double *p,
}
-int get_scalable(Reflection *refl)
+int get_scalable(const Reflection *refl)
{
return refl->data.scalable;
}
-int get_redundancy(Reflection *refl)
+int get_redundancy(const Reflection *refl)
{
return refl->data.redundancy;
}
-double get_sum_squared_dev(Reflection *refl)
+double get_sum_squared_dev(const Reflection *refl)
{
return refl->data.sum_squared_dev;
}
-double get_esd_intensity(Reflection *refl)
+double get_esd_intensity(const Reflection *refl)
{
return refl->data.esd_i;
}
-double get_phase(Reflection *refl)
+double get_phase(const Reflection *refl)
{
return refl->data.phase;
}
@@ -261,7 +262,7 @@ double get_phase(Reflection *refl)
/********************************** Setters ***********************************/
-void copy_data(Reflection *to, Reflection *from)
+void copy_data(Reflection *to, const Reflection *from)
{
memcpy(&to->data, &from->data, sizeof(struct _refldata));
}
diff --git a/src/reflist.h b/src/reflist.h
index 5b647102..462a19af 100644
--- a/src/reflist.h
+++ b/src/reflist.h
@@ -28,26 +28,26 @@ extern RefList *reflist_new(void);
extern void reflist_free(RefList *list);
/* Search */
-extern Reflection *find_refl(RefList *list, INDICES);
+extern Reflection *find_refl(const RefList *list, INDICES);
extern Reflection *next_found_refl(Reflection *refl);
/* Get */
-extern double get_excitation_error(Reflection *refl);
-extern void get_detector_pos(Reflection *refl, double *x, double *y);
-extern void get_indices(Reflection *refl,
+extern double get_excitation_error(const Reflection *refl);
+extern void get_detector_pos(const Reflection *refl, double *x, double *y);
+extern void get_indices(const Reflection *refl,
signed int *h, signed int *k, signed int *l);
-extern double get_partiality(Reflection *refl);
-extern double get_intensity(Reflection *refl);
-extern void get_partial(Reflection *refl, double *r1, double *r2, double *p,
+extern double get_partiality(const Reflection *refl);
+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(Reflection *refl);
-extern int get_redundancy(Reflection *refl);
-extern double get_sum_squared_dev(Reflection *refl);
-extern double get_esd_intensity(Reflection *refl);
-extern double get_phase(Reflection *refl);
+extern int get_scalable(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);
+extern double get_phase(const Reflection *refl);
/* Set */
-extern void copy_data(Reflection *to, Reflection *from);
+extern void copy_data(Reflection *to, const Reflection *from);
extern void set_detector_pos(Reflection *refl, double exerr,
double x, double y);
extern void set_partial(Reflection *refl, double r1, double r2, double p,