aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--doc/man/partialator.111
-rw-r--r--libcrystfel/src/geometry.c11
-rw-r--r--libcrystfel/src/stream.c6
-rw-r--r--libcrystfel/src/xds.c13
-rw-r--r--src/partialator.c112
-rw-r--r--src/post-refinement.c120
-rw-r--r--src/post-refinement.h8
7 files changed, 162 insertions, 119 deletions
diff --git a/doc/man/partialator.1 b/doc/man/partialator.1
index 5c3c24bb..272c97a9 100644
--- a/doc/man/partialator.1
+++ b/doc/man/partialator.1
@@ -62,7 +62,12 @@ Run \fIn\fR cycles of scaling and post refinement.
.PD 0
.IP \fB--no-scale\fR
.PD
-Fix all the overall scaling factors to be unity.
+Disable the scaling part of the refinement calculation.
+
+.PD 0
+.IP \fB--no-pr\fR
+.PD
+Disable the orientation/physics model part of the refinement calculation.
.PD 0
.IP "\fB-m\fR \fImodel\fR"
@@ -149,6 +154,10 @@ partialator -i \fImy.stream \fR-o \fImy.hkl\fR -y \fImypointgroup \fB--model=uni
.PD
partialator -i \fImy.stream \fR-o \fImy.hkl\fR -y \fImypointgroup \fB--model=scsphere --iterations=0\fR
+.IP "Merging with partialities, with scaling but without post-refinement:"
+.PD
+partialator -i \fImy.stream \fR-o \fImy.hkl\fR -y \fImypointgroup \fB--model=scsphere --iterations=1 --no-pr\fR
+
.IP "Merging with partialities, post-refinement and scaling:"
.PD
partialator -i \fImy.stream \fR-o \fImy.hkl\fR -y \fImypointgroup \fB--model=scsphere --iterations=1\fR
diff --git a/libcrystfel/src/geometry.c b/libcrystfel/src/geometry.c
index 048dc9f7..3325923c 100644
--- a/libcrystfel/src/geometry.c
+++ b/libcrystfel/src/geometry.c
@@ -236,6 +236,17 @@ static Reflection *check_reflection(struct image *image, Crystal *cryst,
/* Calculate partiality */
part = partiality(pmodel, rlow, rhigh, pr);
+ if ( isnan(part) ) {
+ ERROR("Assigning NAN partiality!\n");
+ ERROR("%3i %3i %3i rlow = %e, rhigh = %e\n",
+ h, k, l, rlow, rhigh);
+ ERROR("div + m = %e, R = %e, bw = %e\n", del, pr, image->bw);
+ ERROR("D = %e\n", rlow - rhigh);
+ ERROR("klow = %e khigh = %e\n", klow, khigh);
+ ERROR("distance(..) = %e\n", distance(cet, cez, tl, zl));
+ abort();
+ }
+
if ( updateme == NULL ) {
refl = reflection_new(h, k, l);
} else {
diff --git a/libcrystfel/src/stream.c b/libcrystfel/src/stream.c
index 808149e3..313e25e4 100644
--- a/libcrystfel/src/stream.c
+++ b/libcrystfel/src/stream.c
@@ -1151,11 +1151,11 @@ int read_chunk_2(Stream *st, struct image *image, StreamReadFlags srf)
image->div = div;
}
- if ( sscanf(line, "beam_bandwidth = %f %%", &bw) == 1 ) {
- image->bw = bw/100.0;
+ if ( sscanf(line, "beam_bandwidth = %f", &bw) == 1 ) {
+ image->bw = bw;
}
- if ( sscanf(line, "num_peaks = %lld %%", &num_peaks) == 1 ) {
+ if ( sscanf(line, "num_peaks = %lld", &num_peaks) == 1 ) {
image->num_peaks = num_peaks;
}
diff --git a/libcrystfel/src/xds.c b/libcrystfel/src/xds.c
index f5985e73..36898949 100644
--- a/libcrystfel/src/xds.c
+++ b/libcrystfel/src/xds.c
@@ -619,6 +619,19 @@ IndexingPrivate *xds_prepare(IndexingMethod *indm, UnitCell *cell,
struct detector *det, float *ltl)
{
struct xds_private *xp;
+ int need_cell = 0;
+
+ /* Check if cell parameters are needed/provided */
+ if ( *indm & INDEXING_CHECK_CELL_COMBINATIONS ) need_cell = 1;
+ if ( *indm & INDEXING_CHECK_CELL_AXES ) need_cell = 1;
+ if ( *indm & INDEXING_USE_CELL_PARAMETERS ) need_cell = 1;
+ if ( need_cell && !cell_has_parameters(cell) ) {
+ ERROR("Altering your XDS flags because cell parameters were"
+ " not provided.\n");
+ *indm &= ~INDEXING_CHECK_CELL_COMBINATIONS;
+ *indm &= ~INDEXING_CHECK_CELL_AXES;
+ *indm &= ~INDEXING_USE_CELL_PARAMETERS;
+ }
/* Either cell,latt and cell provided, or nocell-nolatt and no cell
* - complain about anything else. Could figure this out automatically,
diff --git a/src/partialator.c b/src/partialator.c
index 46353c25..2ad20094 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -74,7 +74,8 @@ static void show_help(const char *s)
" --start-after=<n> Skip <n> crystals at the start of the stream.\n"
" --stop-after=<n> Stop after merging <n> crystals.\n"
" -n, --iterations=<n> Run <n> cycles of scaling and post-refinement.\n"
-" --no-scale Fix all the scaling factors at unity.\n"
+" --no-scale Disable scale factor (G, B) refinement.\n"
+" --no-pr Disable orientation/physics refinement.\n"
" -m, --model=<model> Specify partiality model.\n"
" --min-measurements=<n> Minimum number of measurements to require.\n"
" --no-polarisation Disable polarisation correction.\n"
@@ -84,111 +85,6 @@ static void show_help(const char *s)
}
-struct refine_args
-{
- RefList *full;
- Crystal *crystal;
- PartialityModel pmodel;
- int no_scale;
- struct prdata prdata;
-};
-
-
-struct queue_args
-{
- int n_started;
- int n_done;
- Crystal **crystals;
- int n_crystals;
- struct refine_args task_defaults;
- double initial_residual;
- double initial_free_residual;
- double final_residual;
- double final_free_residual;
-};
-
-
-static void refine_image(void *task, int id)
-{
- struct refine_args *pargs = task;
- Crystal *cr = pargs->crystal;
-
- pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel,
- pargs->no_scale);
-}
-
-
-static void *get_image(void *vqargs)
-{
- struct refine_args *task;
- struct queue_args *qargs = vqargs;
-
- task = malloc(sizeof(struct refine_args));
- memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
-
- task->crystal = qargs->crystals[qargs->n_started];
-
- qargs->n_started++;
-
- return task;
-}
-
-
-static void done_image(void *vqargs, void *task)
-{
- struct queue_args *qargs = vqargs;
- struct refine_args *pargs = task;
-
- qargs->n_done++;
- qargs->initial_residual += pargs->prdata.initial_residual;
- qargs->initial_free_residual += pargs->prdata.initial_free_residual;
- qargs->final_residual += pargs->prdata.final_residual;
- qargs->final_free_residual += pargs->prdata.final_free_residual;
-
- progress_bar(qargs->n_done, qargs->n_crystals, "Refining");
- free(task);
-}
-
-
-static void refine_all(Crystal **crystals, int n_crystals,
- RefList *full, int nthreads, PartialityModel pmodel,
- int no_scale,
- double *initial_residual, double *initial_free_residual,
- double *final_residual, double *final_free_residual)
-{
- struct refine_args task_defaults;
- struct queue_args qargs;
-
- task_defaults.full = full;
- task_defaults.crystal = NULL;
- task_defaults.pmodel = pmodel;
- task_defaults.prdata.refined = 0;
- task_defaults.prdata.n_filtered = 0;
- task_defaults.no_scale = no_scale;
-
- qargs.task_defaults = task_defaults;
- qargs.n_started = 0;
- qargs.n_done = 0;
- qargs.n_crystals = n_crystals;
- qargs.crystals = crystals;
- qargs.initial_residual = 0.0;
- qargs.initial_free_residual = 0.0;
- qargs.final_residual = 0.0;
- qargs.final_free_residual = 0.0;
-
- /* Don't have threads which are doing nothing */
- if ( n_crystals < nthreads ) nthreads = n_crystals;
-
- run_threads(nthreads, refine_image, get_image, done_image,
- &qargs, n_crystals, 0, 0, 0);
-
- *initial_residual = qargs.initial_residual;
- *initial_free_residual = qargs.initial_free_residual;
- *final_residual = qargs.final_residual;
- *final_free_residual = qargs.final_free_residual;
-}
-
-
static void display_progress(int n_images, int n_crystals)
{
if ( !isatty(STDERR_FILENO) ) return;
@@ -475,6 +371,7 @@ int main(int argc, char *argv[])
int n_crystals_seen = 0;
char cmdline[1024];
int no_scale = 0;
+ int no_pr = 0;
Stream *st;
Crystal **crystals;
char *pmodel_str = NULL;
@@ -510,6 +407,7 @@ int main(int argc, char *argv[])
{"res-push", 1, NULL, 5}, /* compat */
{"no-scale", 0, &no_scale, 1},
+ {"no-pr", 0, &no_pr, 1},
{"no-polarisation", 0, &polarisation, 0},
{"no-polarization", 0, &polarisation, 0},
{"polarisation", 0, &polarisation, 1}, /* compat */
@@ -831,7 +729,7 @@ int main(int argc, char *argv[])
/* Refine all crystals to get the best fit */
refine_all(crystals, n_crystals, full, nthreads, pmodel,
- no_scale, &init_dev, &init_free_dev,
+ no_scale, no_pr, &init_dev, &init_free_dev,
&final_dev, &final_free_dev);
STATUS("Overall residual: initial = %e, final = %e\n",
diff --git a/src/post-refinement.c b/src/post-refinement.c
index 139b8b9a..d62d9323 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -326,8 +326,8 @@ static void apply_shift(Crystal *cr, int k, double shift)
/* Perform one cycle of post refinement on 'image' against 'full' */
static double pr_iterate(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int no_scale, int *n_filtered,
- int verbose)
+ PartialityModel pmodel, int no_scale, int no_pr,
+ int *n_filtered, int verbose)
{
gsl_matrix *M;
gsl_vector *v;
@@ -346,7 +346,7 @@ static double pr_iterate(Crystal *cr, const RefList *full,
/* If partiality model is anything other than "unity", refine all the
* geometrical parameters */
- if ( pmodel != PMODEL_UNITY ) {
+ if ( (pmodel != PMODEL_UNITY) && !no_pr ) {
rv[num_params++] = GPARAM_ASX;
rv[num_params++] = GPARAM_ASY;
rv[num_params++] = GPARAM_ASZ;
@@ -546,8 +546,8 @@ static double residual(Crystal *cr, const RefList *full, int verbose, int free)
}
-struct prdata pr_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int no_scale)
+static struct prdata pr_refine(Crystal *cr, const RefList *full,
+ PartialityModel pmodel, int no_scale, int no_pr)
{
int i;
int verbose = 0;
@@ -584,7 +584,8 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
cell_get_reciprocal(crystal_get_cell(cr), &asx, &asy, &asz,
&bsx, &bsy, &bsz, &csx, &csy, &csz);
- pr_iterate(cr, full, pmodel, no_scale, &prdata.n_filtered, 0);
+ pr_iterate(cr, full, pmodel, no_scale, no_pr,
+ &prdata.n_filtered, verbose);
update_partialities(cr, pmodel);
@@ -614,3 +615,110 @@ struct prdata pr_refine(Crystal *cr, const RefList *full,
return prdata;
}
+
+
+struct refine_args
+{
+ RefList *full;
+ Crystal *crystal;
+ PartialityModel pmodel;
+ int no_scale;
+ int no_pr;
+ struct prdata prdata;
+};
+
+
+struct queue_args
+{
+ int n_started;
+ int n_done;
+ Crystal **crystals;
+ int n_crystals;
+ struct refine_args task_defaults;
+ double initial_residual;
+ double initial_free_residual;
+ double final_residual;
+ double final_free_residual;
+};
+
+
+static void refine_image(void *task, int id)
+{
+ struct refine_args *pargs = task;
+ Crystal *cr = pargs->crystal;
+
+ pargs->prdata = pr_refine(cr, pargs->full, pargs->pmodel,
+ pargs->no_scale, pargs->no_pr);
+}
+
+
+static void *get_image(void *vqargs)
+{
+ struct refine_args *task;
+ struct queue_args *qargs = vqargs;
+
+ task = malloc(sizeof(struct refine_args));
+ memcpy(task, &qargs->task_defaults, sizeof(struct refine_args));
+
+ task->crystal = qargs->crystals[qargs->n_started];
+
+ qargs->n_started++;
+
+ return task;
+}
+
+
+static void done_image(void *vqargs, void *task)
+{
+ struct queue_args *qargs = vqargs;
+ struct refine_args *pargs = task;
+
+ qargs->n_done++;
+ qargs->initial_residual += pargs->prdata.initial_residual;
+ qargs->initial_free_residual += pargs->prdata.initial_free_residual;
+ qargs->final_residual += pargs->prdata.final_residual;
+ qargs->final_free_residual += pargs->prdata.final_free_residual;
+
+ progress_bar(qargs->n_done, qargs->n_crystals, "Refining");
+ free(task);
+}
+
+
+void refine_all(Crystal **crystals, int n_crystals,
+ RefList *full, int nthreads, PartialityModel pmodel,
+ int no_scale, int no_pr,
+ double *initial_residual, double *initial_free_residual,
+ double *final_residual, double *final_free_residual)
+{
+ struct refine_args task_defaults;
+ struct queue_args qargs;
+
+ task_defaults.full = full;
+ task_defaults.crystal = NULL;
+ task_defaults.pmodel = pmodel;
+ task_defaults.prdata.refined = 0;
+ task_defaults.prdata.n_filtered = 0;
+ task_defaults.no_scale = no_scale;
+ task_defaults.no_pr = no_pr;
+
+ qargs.task_defaults = task_defaults;
+ qargs.n_started = 0;
+ qargs.n_done = 0;
+ qargs.n_crystals = n_crystals;
+ qargs.crystals = crystals;
+ qargs.initial_residual = 0.0;
+ qargs.initial_free_residual = 0.0;
+ qargs.final_residual = 0.0;
+ qargs.final_free_residual = 0.0;
+
+ /* Don't have threads which are doing nothing */
+ if ( n_crystals < nthreads ) nthreads = n_crystals;
+
+ run_threads(nthreads, refine_image, get_image, done_image,
+ &qargs, n_crystals, 0, 0, 0);
+
+ *initial_residual = qargs.initial_residual;
+ *initial_free_residual = qargs.initial_free_residual;
+ *final_residual = qargs.final_residual;
+ *final_free_residual = qargs.final_free_residual;
+}
diff --git a/src/post-refinement.h b/src/post-refinement.h
index fd2d771b..b59ae24f 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -54,8 +54,12 @@ struct prdata
};
-extern struct prdata pr_refine(Crystal *cr, const RefList *full,
- PartialityModel pmodel, int no_scale);
+
+extern void refine_all(Crystal **crystals, int n_crystals,
+ RefList *full, int nthreads, PartialityModel pmodel,
+ int no_scale, int no_pr,
+ double *initial_residual, double *initial_free_residual,
+ double *final_residual, double *final_free_residual);
/* Exported so it can be poked by tests/pr_p_gradient_check */
extern double gradient(Crystal *cr, int k, Reflection *refl,
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-02 07:20:17 +0000