aboutsummaryrefslogtreecommitdiff
path: root/src/post-refinement.c
diff options
context:
space:
mode:
Diffstat (limited to 'src/post-refinement.c')
-rw-r--r--src/post-refinement.c185
1 files changed, 92 insertions, 93 deletions
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 70b14a0f..4a02ed08 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -56,14 +56,16 @@ struct rf_alteration
struct rf_priv
{
- Crystal *cr; /**< Original crystal (before any refinement) */
+ Crystal *cr; /**< Original crystal (before any refinement) */
+ struct image *image; /**< Original image (before any refinement) */
const RefList *full;
int serial;
int scaleflags;
PartialityModel pmodel;
- Crystal *cr_tgt; /**< Crystal to use for testing modifications */
- struct image image_tgt; /**< Image structure to go with cr_tgt */
+ Crystal *cr_tgt; /**< Crystal to use for testing modifications */
+ struct image image_tgt; /**< Image structure to go with cr_tgt */
+ RefList *refls; /**< The reflections to use */
};
@@ -151,32 +153,30 @@ static void rotate_cell_xy(UnitCell *source, UnitCell *tgt,
static void apply_parameters(Crystal *cr_orig, Crystal *cr_tgt,
+ struct image *im_orig, struct image *im_tgt,
struct rf_alteration alter)
{
- double R, lambda;
+ double R;
struct gaussian g;
- struct image *image;
- image = crystal_get_image(cr_tgt);
R = crystal_get_profile_radius(cr_orig);
- lambda = crystal_get_image_const(cr_orig)->lambda;
rotate_cell_xy(crystal_get_cell(cr_orig), crystal_get_cell(cr_tgt),
alter.rot_x, alter.rot_y);
crystal_set_profile_radius(cr_tgt, R+alter.delta_R);
- image->lambda = lambda+alter.delta_wave;
+ im_tgt->lambda = im_orig->lambda+alter.delta_wave;
- g.kcen = 1.0/image->lambda;
- g.sigma = image->bw/image->lambda;
+ g.kcen = 1.0/im_tgt->lambda;
+ g.sigma = im_tgt->bw/im_tgt->lambda;
g.area = 1;
- spectrum_set_gaussians(image->spectrum, &g, 1);
+ spectrum_set_gaussians(im_tgt->spectrum, &g, 1);
}
static double calc_residual(struct rf_priv *pv, struct rf_alteration alter,
int free)
{
- apply_parameters(pv->cr, pv->cr_tgt, alter);
+ apply_parameters(pv->cr, pv->cr_tgt, pv->image, &pv->image_tgt, alter);
if ( fabs(crystal_get_profile_radius(pv->cr_tgt)) > 5e9 ) {
STATUS("radius > 5e9\n");
@@ -189,16 +189,16 @@ static double calc_residual(struct rf_priv *pv, struct rf_alteration alter,
return NAN;
}
- if ( crystal_get_image(pv->cr_tgt)->lambda <= 0.0 ) {
+ if ( pv->image_tgt.lambda <= 0.0 ) {
STATUS("lambda < 0\n");
return NAN;
}
- update_predictions(pv->cr_tgt);
- calculate_partialities(pv->cr_tgt, pv->pmodel);
+ update_predictions(pv->refls, pv->cr_tgt, &pv->image_tgt);
+ calculate_partialities(pv->refls, pv->cr_tgt, &pv->image_tgt, pv->pmodel);
- return residual(pv->cr_tgt, pv->full, free, NULL, NULL);
+ return residual(pv->refls, pv->cr_tgt, pv->full, free, NULL, NULL);
}
@@ -264,64 +264,58 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx)
}
-static void try_reindex(Crystal *crin, const RefList *full,
- SymOpList *sym, SymOpList *amb, int scaleflags,
- PartialityModel pmodel)
+static void try_reindex(RefList **prefls, Crystal *crin, struct image *image,
+ const RefList *full, SymOpList *sym, SymOpList *amb,
+ int scaleflags, PartialityModel pmodel)
{
- Crystal *cr;
double residual_original;
int idx, n;
if ( sym == NULL || amb == NULL ) return;
- if ( scale_one_crystal(crin, full, scaleflags) ) return;
- residual_original = residual(crin, full, 0, NULL, NULL);
-
- cr = crystal_copy(crin);
+ if ( scale_one_crystal(*prefls, crin, full, scaleflags) ) return;
+ residual_original = residual(*prefls, crin, full, 0, NULL, NULL);
n = num_equivs(amb, NULL);
for ( idx=0; idx<n; idx++ ) {
RefList *list;
- UnitCell *cell;
- double residual_flipped;
-
- cell = cell_new_from_cell(crystal_get_cell(crin));
- if ( cell == NULL ) return;
- reindex_cell(cell, amb, idx);
- crystal_set_cell(cr, cell);
-
- list = reindex_reflections(crystal_get_reflections(crin),
- amb, sym, idx);
- crystal_set_reflections(cr, list);
-
- update_predictions(cr);
- calculate_partialities(cr, pmodel);
-
- if ( scale_one_crystal(cr, full, scaleflags) ) return;
- residual_flipped = residual(cr, full, 0, NULL, NULL);
-
- if ( residual_flipped < residual_original ) {
- crystal_set_cell(crin, cell);
- crystal_set_reflections(crin, list);
- residual_original = residual_flipped;
- } else {
- cell_free(crystal_get_cell(cr));
- reflist_free(crystal_get_reflections(cr));
+ Crystal *cr;
+
+ cr = crystal_copy(crin);
+
+ reindex_cell(crystal_get_cell(cr), amb, idx);
+ list = reindex_reflections(*prefls, amb, sym, idx);
+
+ update_predictions(list, cr, image);
+ calculate_partialities(list, cr, image, pmodel);
+ if ( scale_one_crystal(list, cr, full, scaleflags) == 0 ) {
+
+ double residual_flipped;
+ residual_flipped = residual(list, cr, full, 0, NULL, NULL);
+
+ if ( residual_flipped < residual_original ) {
+ crystal_set_cell(crin, crystal_relinquish_cell(cr));
+ reflist_free(*prefls);
+ *prefls = list;
+ residual_original = residual_flipped;
+ } else {
+ reflist_free(list);
+ }
}
+
+ crystal_free(cr);
}
- crystal_free(cr);
}
-void write_test_logs(Crystal *crystal, const RefList *full,
+void write_test_logs(Crystal *crystal, struct image *image, const RefList *full,
signed int cycle, int serial,
const char *log_folder)
{
FILE *fh;
- struct image *image = crystal_get_image(crystal);
char tmp[256];
char ins[16];
@@ -367,7 +361,7 @@ void write_test_logs(Crystal *crystal, const RefList *full,
}
-void write_specgraph(Crystal *crystal, const RefList *full,
+void write_specgraph(RefList *list, Crystal *crystal, struct image *image, const RefList *full,
signed int cycle, int serial,
const char *log_folder)
{
@@ -378,7 +372,6 @@ void write_specgraph(Crystal *crystal, const RefList *full,
double G = crystal_get_osf(crystal);
double B = crystal_get_Bfac(crystal);
UnitCell *cell;
- struct image *image = crystal_get_image(crystal);
char ins[16];
snprintf(tmp, 256, "%s/specgraph-crystal%i.dat", log_folder, serial);
@@ -408,7 +401,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
ins[1] = '\0';
}
- for ( refl = first_refl(crystal_get_reflections(crystal), &iter);
+ for ( refl = first_refl(list, &iter);
refl != NULL;
refl = next_refl(refl, iter) )
{
@@ -443,7 +436,8 @@ void write_specgraph(Crystal *crystal, const RefList *full,
}
-static void write_angle_grid(Crystal *cr, const RefList *full,
+static void write_angle_grid(RefList *list_in, Crystal *cr, struct image *image,
+ const RefList *full,
signed int cycle, int serial, int scaleflags,
PartialityModel pmodel,
const char *log_folder)
@@ -452,7 +446,6 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
char fn[64];
char ins[16];
struct rf_priv priv;
- RefList *list;
UnitCell *cell;
Spectrum *spectrum;
@@ -462,12 +455,11 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
priv.scaleflags = scaleflags;
priv.pmodel = pmodel;
priv.cr_tgt = crystal_copy(cr);
- priv.image_tgt = *crystal_get_image(cr);
+ priv.image = image;
+ priv.image_tgt = *image;
spectrum = spectrum_new();
priv.image_tgt.spectrum = spectrum;
- crystal_set_image(priv.cr_tgt, &priv.image_tgt);
- list = copy_reflist(crystal_get_reflections(cr));
- crystal_set_reflections(priv.cr_tgt, list);
+ priv.refls = copy_reflist(list_in);
cell = cell_new_from_cell(crystal_get_cell(cr));
crystal_set_cell(priv.cr_tgt, cell);
@@ -504,13 +496,14 @@ static void write_angle_grid(Crystal *cr, const RefList *full,
fclose(fh);
}
- reflist_free(crystal_get_reflections(priv.cr_tgt));
+ reflist_free(priv.refls);
crystal_free(priv.cr_tgt);
spectrum_free(spectrum);
}
-static void write_radius_grid(Crystal *cr, const RefList *full,
+static void write_radius_grid(RefList *list_in, Crystal *cr, struct image *image,
+ const RefList *full,
signed int cycle, int serial, int scaleflags,
PartialityModel pmodel,
const char *log_folder)
@@ -519,7 +512,6 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
char fn[64];
char ins[16];
struct rf_priv priv;
- RefList *list;
UnitCell *cell;
Spectrum *spectrum;
@@ -529,12 +521,11 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
priv.scaleflags = scaleflags;
priv.pmodel = pmodel;
priv.cr_tgt = crystal_copy(cr);
- priv.image_tgt = *crystal_get_image(cr);
+ priv.image = image;
+ priv.image_tgt = *image;
spectrum = spectrum_new();
priv.image_tgt.spectrum = spectrum;
- crystal_set_image(priv.cr_tgt, &priv.image_tgt);
- list = copy_reflist(crystal_get_reflections(cr));
- crystal_set_reflections(priv.cr_tgt, list);
+ priv.refls = copy_reflist(list_in);
cell = cell_new_from_cell(crystal_get_cell(cr));
crystal_set_cell(priv.cr_tgt, cell);
@@ -571,19 +562,20 @@ static void write_radius_grid(Crystal *cr, const RefList *full,
fclose(fh);
}
- reflist_free(crystal_get_reflections(priv.cr_tgt));
+ reflist_free(priv.refls);
crystal_free(priv.cr_tgt);
spectrum_free(spectrum);
}
-void write_gridscan(Crystal *cr, const RefList *full,
+void write_gridscan(RefList *list,
+ Crystal *cr, struct image *image, const RefList *full,
signed int cycle, int serial, int scaleflags,
PartialityModel pmodel,
const char *log_folder)
{
- write_angle_grid(cr, full, cycle, serial, scaleflags, pmodel, log_folder);
- write_radius_grid(cr, full, cycle, serial, scaleflags, pmodel, log_folder);
+ write_angle_grid(list, cr, image, full, cycle, serial, scaleflags, pmodel, log_folder);
+ write_radius_grid(list, cr, image, full, cycle, serial, scaleflags, pmodel, log_folder);
}
@@ -631,7 +623,7 @@ static int refine_loop(struct rf_alteration *cur, struct rf_alteration *dirns,
*total_iter, lowest_fom, freefom,
cur->rot_x, cur->rot_y,
crystal_get_profile_radius(priv->cr)+cur->delta_R,
- crystal_get_image_const(priv->cr)->lambda+cur->delta_wave);
+ priv->image->lambda+cur->delta_wave);
}
} while ( (best != 0) && (n_iter < 10) );
@@ -649,7 +641,8 @@ static void zero_alter(struct rf_alteration *alter)
}
-static void do_pr_refine(Crystal *cr, const RefList *full,
+static void do_pr_refine(RefList **plist_in, Crystal *cr, struct image *image,
+ const RefList *full,
PartialityModel pmodel, int serial,
int cycle, int write_logs,
SymOpList *sym, SymOpList *amb, int scaleflags,
@@ -659,12 +652,11 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
struct rf_alteration alter;
int n_iter = 0;
double fom, freefom;
- RefList *list;
FILE *fh = NULL;
UnitCell *cell;
Spectrum *spectrum;
- try_reindex(cr, full, sym, amb, scaleflags, pmodel);
+ try_reindex(plist_in, cr, image, full, sym, amb, scaleflags, pmodel);
zero_alter(&alter);
@@ -674,12 +666,11 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
priv.scaleflags = scaleflags;
priv.pmodel = pmodel;
priv.cr_tgt = crystal_copy(cr);
- priv.image_tgt = *crystal_get_image(cr);
+ priv.image = image;
+ priv.image_tgt = *image;
spectrum = spectrum_new();
priv.image_tgt.spectrum = spectrum;
- crystal_set_image(priv.cr_tgt, &priv.image_tgt);
- list = copy_reflist(crystal_get_reflections(cr));
- crystal_set_reflections(priv.cr_tgt, list);
+ priv.refls = copy_reflist(*plist_in);
cell = cell_new_from_cell(crystal_get_cell(cr));
crystal_set_cell(priv.cr_tgt, cell);
@@ -700,7 +691,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
n_iter, fom, freefom,
alter.rot_x, alter.rot_y,
crystal_get_profile_radius(cr)+alter.delta_R,
- crystal_get_image(cr)->lambda+alter.delta_wave);
+ image->lambda+alter.delta_wave);
}
}
@@ -724,15 +715,15 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
refine_loop(&alter, dirns, 2, &priv, &n_iter, fh);
/* Apply the final shifts */
- apply_parameters(cr, cr, alter);
- update_predictions(cr);
- calculate_partialities(cr, pmodel);
+ apply_parameters(cr, cr, image, image, alter);
+ update_predictions(*plist_in, cr, image);
+ calculate_partialities(*plist_in, cr, image, pmodel);
if ( write_logs ) {
- write_gridscan(cr, full, cycle, serial, scaleflags,
+ write_gridscan(*plist_in, cr, image, full, cycle, serial, scaleflags,
pmodel, log_folder);
- write_specgraph(cr, full, cycle, serial, log_folder);
- write_test_logs(cr, full, cycle, serial, log_folder);
+ write_specgraph(*plist_in, cr, image, full, cycle, serial, log_folder);
+ write_test_logs(cr, image, full, cycle, serial, log_folder);
}
if ( crystal_get_profile_radius(cr) > 5e9 ) {
@@ -743,7 +734,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
fclose(fh);
}
- reflist_free(crystal_get_reflections(priv.cr_tgt));
+ reflist_free(priv.refls);
crystal_free(priv.cr_tgt);
spectrum_free(spectrum);
}
@@ -753,6 +744,8 @@ struct refine_args
{
RefList *full;
Crystal *crystal;
+ RefList **prefls;
+ struct image *image;
PartialityModel pmodel;
int serial;
int cycle;
@@ -768,7 +761,8 @@ struct pr_queue_args
{
int n_started;
int n_done;
- Crystal **crystals;
+ struct crystal_refls *crystals;
+ struct image **images;
int n_crystals;
struct refine_args task_defaults;
};
@@ -777,12 +771,12 @@ struct pr_queue_args
static void refine_image(void *task, int id)
{
struct refine_args *pargs = task;
- Crystal *cr = pargs->crystal;
int write_logs = 0;
write_logs = !pargs->no_logs && (pargs->serial % 20 == 0);
- do_pr_refine(cr, pargs->full, pargs->pmodel,
+ do_pr_refine(pargs->prefls, pargs->crystal, pargs->image,
+ pargs->full, pargs->pmodel,
pargs->serial, pargs->cycle, write_logs,
pargs->sym, pargs->amb, pargs->scaleflags,
pargs->log_folder);
@@ -797,7 +791,9 @@ static void *get_image(void *vqargs)
task = malloc(sizeof(struct refine_args));
memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
- task->crystal = qargs->crystals[qargs->n_started];
+ task->crystal = qargs->crystals[qargs->n_started].cr;
+ task->prefls = &qargs->crystals[qargs->n_started].refls;
+ task->image = qargs->images[qargs->n_started];
task->serial = qargs->n_started;
qargs->n_started++;
@@ -817,7 +813,7 @@ static void done_image(void *vqargs, void *task)
}
-void refine_all(Crystal **crystals, int n_crystals,
+void refine_all(struct crystal_refls *crystals, struct image **images, int n_crystals,
RefList *full, int nthreads, PartialityModel pmodel,
int cycle, int no_logs,
SymOpList *sym, SymOpList *amb, int scaleflags,
@@ -828,6 +824,8 @@ void refine_all(Crystal **crystals, int n_crystals,
task_defaults.full = full;
task_defaults.crystal = NULL;
+ task_defaults.prefls = NULL;
+ task_defaults.image = NULL;
task_defaults.pmodel = pmodel;
task_defaults.cycle = cycle;
task_defaults.no_logs = no_logs;
@@ -842,6 +840,7 @@ void refine_all(Crystal **crystals, int n_crystals,
qargs.n_done = 0;
qargs.n_crystals = n_crystals;
qargs.crystals = crystals;
+ qargs.images = images;
/* Don't have threads which are doing nothing */
if ( n_crystals < nthreads ) nthreads = n_crystals;
generated by cgit v1.2.3 (git 2.25.1) at 2024-05-17 05:45:36 +0000