Skip to content
GitLab
Commit 656abb22 authored by Thomas White's avatar Thomas White
Browse files

partialator: Save reflection list

parent 407a74d5
Hide whitespace changes
Inline Side-by-side
src/hrs-scaling.c
View file @ 656abb22
......@@ -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;
}
src/hrs-scaling.h
View file @ 656abb22
......@@ -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);
......
src/partialator.c
View file @ 656abb22
......@@ -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);
......
src/post-refinement.c
View file @ 656abb22
......@@ -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) );
......
src/post-refinement.h
View file @ 656abb22
......@@ -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);
......
src/reflist.c
View file @ 656abb22
......@@ -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));
}
......
src/reflist.h
View file @ 656abb22
......@@ -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,
......
Supports Markdown
0% or .
You are about to add 0 people to the discussion. Proceed with caution.
Finish editing this message first!
Please register or to comment

imprint | privacy | Conditions of Use