aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/partialator.c35
-rw-r--r--src/post-refinement.c55
-rw-r--r--src/post-refinement.h4
-rw-r--r--src/scaling.c27
-rw-r--r--src/scaling.h18
-rw-r--r--tests/linear_scale_check.c10
6 files changed, 99 insertions, 50 deletions
diff --git a/src/partialator.c b/src/partialator.c
index 0f3bd104..e775caef 100644
--- a/src/partialator.c
+++ b/src/partialator.c
@@ -317,6 +317,7 @@ static void show_help(const char *s)
" --stop-after=<n> Stop after merging <n> crystals.\n"
" -n, --iterations=<n> Run <n> cycles of scaling and post-refinement.\n"
" --no-scale Disable scale factor (G, B) refinement.\n"
+" --no-Bscale Disable B factor scaling.\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"
@@ -717,6 +718,7 @@ static void all_residuals(Crystal **crystals, int n_crystals, RefList *full,
if ( crystal_get_user_flag(crystals[i]) ) continue;
+ /* Scaling should have been done right before calling this */
r = residual(crystals[i], full, 0, NULL, NULL);
free_r = residual(crystals[i], full, 1, NULL, NULL);
log_r = log_residual(crystals[i], full, 0, NULL, NULL);
@@ -754,6 +756,7 @@ struct log_qargs
Crystal **crystals;
int n_crystals;
RefList *full;
+ int scaleflags;
};
@@ -761,6 +764,7 @@ struct log_args
{
Crystal *cr;
RefList *full;
+ int scaleflags;
int iter;
int cnum;
};
@@ -780,6 +784,7 @@ static void *get_log_task(void *vp)
task->full = qargs->full;
task->iter = qargs->iter;
task->cnum = qargs->next;
+ task->scaleflags = qargs->scaleflags;
qargs->next += 20;
return task;
@@ -790,7 +795,8 @@ static void write_logs(void *vp, int cookie)
{
struct log_args *args = vp;
write_specgraph(args->cr, args->full, args->iter, args->cnum);
- write_gridscan(args->cr, args->full, args->iter, args->cnum);
+ write_gridscan(args->cr, args->full, args->iter, args->cnum,
+ args->scaleflags);
write_test_logs(args->cr, args->full, args->iter, args->cnum);
}
@@ -803,16 +809,17 @@ static void done_log(void *qargs, void *vp)
static void write_logs_parallel(Crystal **crystals, int n_crystals,
- RefList *full, int iter, int n_threads)
+ RefList *full, int iter, int n_threads,
+ int scaleflags)
{
struct log_qargs qargs;
-
qargs.iter = iter;
qargs.next = 0;
qargs.full = full;
qargs.crystals = crystals;
qargs.n_crystals = n_crystals;
+ qargs.scaleflags = scaleflags;
STATUS("Writing logs...\n");
run_threads(n_threads, write_logs, get_log_task, done_log, &qargs,
@@ -839,6 +846,7 @@ int main(int argc, char *argv[])
int n_crystals_seen = 0;
char cmdline[1024];
int no_scale = 0;
+ int no_Bscale = 0;
int no_pr = 0;
Stream *st;
Crystal **crystals;
@@ -867,6 +875,7 @@ int main(int argc, char *argv[])
double force_bandwidth = -1.0;
double force_radius = -1.0;
char *audit_info;
+ int scaleflags = 0;
/* Long options */
const struct option longopts[] = {
@@ -895,6 +904,7 @@ int main(int argc, char *argv[])
{"force-radius", 1, NULL, 11},
{"no-scale", 0, &no_scale, 1},
+ {"no-Bscale", 0, &no_Bscale, 1},
{"no-pr", 0, &no_pr, 1},
{"no-polarisation", 0, &polarisation, 0},
{"no-polarization", 0, &polarisation, 0}, /* compat */
@@ -1146,6 +1156,10 @@ int main(int argc, char *argv[])
"partialities (--model=unity).\n");
}
+ if ( no_Bscale ) {
+ scaleflags |= SCALE_NO_B;
+ }
+
if ( !no_logs ) {
int r = mkdir("pr-logs", S_IRWXU | S_IRWXG | S_IROTH | S_IXOTH);
if ( r ) {
@@ -1327,7 +1341,7 @@ int main(int argc, char *argv[])
/* Create reference data set if we don't already have one */
if ( reference == NULL ) {
STATUS("Initial scaling...\n");
- scale_all(crystals, n_crystals, nthreads);
+ scale_all(crystals, n_crystals, nthreads, scaleflags);
full = merge_intensities(crystals, n_crystals, nthreads,
min_measurements, push_res, 1);
} else {
@@ -1339,7 +1353,8 @@ int main(int argc, char *argv[])
show_all_residuals(crystals, n_crystals, full);
if ( !no_pr && !no_logs ) {
write_pgraph(full, crystals, n_crystals, 0, "");
- write_logs_parallel(crystals, n_crystals, full, 0, nthreads);
+ write_logs_parallel(crystals, n_crystals, full, 0, nthreads,
+ scaleflags);
}
/* Iterate */
@@ -1349,14 +1364,15 @@ int main(int argc, char *argv[])
if ( !no_pr ) {
refine_all(crystals, n_crystals, full, nthreads, pmodel,
- 0, i+1, no_logs, sym, amb);
+ 0, i+1, no_logs, sym, amb, scaleflags);
}
/* Create new reference if needed */
if ( reference == NULL ) {
reflist_free(full);
if ( !no_scale ) {
- scale_all(crystals, n_crystals, nthreads);
+ scale_all(crystals, n_crystals, nthreads,
+ scaleflags);
}
full = merge_intensities(crystals, n_crystals, nthreads,
min_measurements,
@@ -1403,7 +1419,7 @@ int main(int argc, char *argv[])
if ( reference == NULL ) {
reflist_free(full);
if ( !no_scale ) {
- scale_all(crystals, n_crystals, nthreads);
+ scale_all(crystals, n_crystals, nthreads, scaleflags);
}
full = merge_intensities(crystals, n_crystals, nthreads,
min_measurements,
@@ -1417,7 +1433,8 @@ int main(int argc, char *argv[])
show_all_residuals(crystals, n_crystals, full);
if ( !no_pr && !no_logs ) {
write_pgraph(full, crystals, n_crystals, -1, "");
- write_logs_parallel(crystals, n_crystals, full, -1, nthreads);
+ write_logs_parallel(crystals, n_crystals, full, -1, nthreads,
+ scaleflags);
}
/* Output results */
diff --git a/src/post-refinement.c b/src/post-refinement.c
index f260da7f..293fe332 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -95,10 +95,6 @@ double residual(Crystal *cr, const RefList *full, int free,
double B = crystal_get_Bfac(cr);
UnitCell *cell = crystal_get_cell(cr);
- if ( linear_scale(full, crystal_get_reflections(cr), &G, complain) ) {
- return GSL_NAN;
- }
-
for ( refl = first_refl(crystal_get_reflections(cr), &iter);
refl != NULL;
refl = next_refl(refl, iter) )
@@ -232,6 +228,7 @@ struct rf_priv {
int verbose;
int serial;
gsl_vector *initial;
+ int scaleflags;
/* For freeing later */
gsl_vector *vals;
@@ -349,6 +346,11 @@ static double residual_f(const gsl_vector *v, void *pp)
update_predictions(cr);
calculate_partialities(cr, PMODEL_XSPHERE);
+ if ( scale_one_crystal(cr, pv->full, pv->scaleflags) ) {
+ crystal_free(cr);
+ if ( pv->verbose ) STATUS("Bad scaling\n");
+ return GSL_NAN;
+ }
res = residual(cr, pv->full, 0, NULL, NULL);
cell_free(crystal_get_cell(cr));
@@ -463,8 +465,8 @@ static void reindex_cell(UnitCell *cell, SymOpList *amb, int idx)
}
-void try_reindex(Crystal *crin, const RefList *full,
- SymOpList *sym, SymOpList *amb)
+static void try_reindex(Crystal *crin, const RefList *full,
+ SymOpList *sym, SymOpList *amb, int scaleflags)
{
RefList *list;
Crystal *cr;
@@ -474,6 +476,7 @@ void try_reindex(Crystal *crin, const RefList *full,
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);
@@ -494,6 +497,7 @@ void try_reindex(Crystal *crin, const RefList *full,
update_predictions(cr);
calculate_partialities(cr, PMODEL_XSPHERE);
+ if ( scale_one_crystal(cr, full, scaleflags) ) return;
residual_flipped = residual(cr, full, 0, NULL, NULL);
if ( residual_flipped < residual_original ) {
@@ -637,6 +641,7 @@ void write_specgraph(Crystal *crystal, const RefList *full,
static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full,
int verbose, int serial,
+ int scaleflags,
struct rf_priv *priv)
{
gsl_multimin_fminimizer *min;
@@ -654,6 +659,7 @@ static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full
priv->full = full;
priv->verbose = verbose;
priv->serial = serial;
+ priv->scaleflags = scaleflags;
priv->f.f = residual_f;
priv->f.n = n_params;
@@ -678,7 +684,7 @@ static gsl_multimin_fminimizer *setup_minimiser(Crystal *cr, const RefList *full
static void write_grid(Crystal *cr, const RefList *full,
- signed int cycle, int serial,
+ signed int cycle, int serial, int scaleflags,
const enum gparam par1, const enum gparam par2,
const char *name)
{
@@ -690,7 +696,7 @@ static void write_grid(Crystal *cr, const RefList *full,
int idx1, idx2;
int i;
- min = setup_minimiser(cr, full, 0, serial, &priv);
+ min = setup_minimiser(cr, full, 0, serial, scaleflags, &priv);
idx1 = 99;
idx2 = 99;
@@ -747,13 +753,13 @@ static void write_grid(Crystal *cr, const RefList *full,
void write_gridscan(Crystal *cr, const RefList *full,
- signed int cycle, int serial)
+ signed int cycle, int serial, int scaleflags)
{
- write_grid(cr, full, cycle, serial,
+ write_grid(cr, full, cycle, serial, scaleflags,
GPARAM_ANG1, GPARAM_ANG2, "ang1-ang2");
- write_grid(cr, full, cycle, serial,
+ write_grid(cr, full, cycle, serial, scaleflags,
GPARAM_ANG1, GPARAM_WAVELENGTH, "ang1-wave");
- write_grid(cr, full, cycle, serial,
+ write_grid(cr, full, cycle, serial, scaleflags,
GPARAM_R, GPARAM_WAVELENGTH, "R-wave");
}
@@ -761,19 +767,21 @@ void write_gridscan(Crystal *cr, const RefList *full,
static void do_pr_refine(Crystal *cr, const RefList *full,
PartialityModel pmodel, int verbose, int serial,
int cycle, int write_logs,
- SymOpList *sym, SymOpList *amb)
+ SymOpList *sym, SymOpList *amb, int scaleflags)
{
gsl_multimin_fminimizer *min;
struct rf_priv priv;
int n_iter = 0;
int status;
- int r;
- double G;
double residual_start, residual_free_start;
FILE *fh = NULL;
- try_reindex(cr, full, sym, amb);
+ try_reindex(cr, full, sym, amb, scaleflags);
+ if ( scale_one_crystal(cr, full, scaleflags) ) {
+ ERROR("Bad scaling at start of refinement.\n");
+ return;
+ }
residual_start = residual(cr, full, 0, NULL, NULL);
residual_free_start = residual(cr, full, 1, NULL, NULL);
@@ -782,7 +790,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
residual_start, residual_free_start);
}
- min = setup_minimiser(cr, full, verbose, serial, &priv);
+ min = setup_minimiser(cr, full, verbose, serial, scaleflags, &priv);
if ( verbose ) {
double res = residual_f(min->x, &priv);
@@ -902,10 +910,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
apply_parameters(min->x, priv.initial, priv.rv, cr);
update_predictions(cr);
calculate_partialities(cr, PMODEL_XSPHERE);
- r = linear_scale(full, crystal_get_reflections(cr), &G, 0);
- if ( r == 0 ) {
- crystal_set_osf(cr, G);
- }
+ scale_one_crystal(cr, full, scaleflags);
if ( verbose ) {
@@ -922,7 +927,7 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
}
if ( write_logs ) {
- write_gridscan(cr, full, cycle, serial);
+ write_gridscan(cr, full, cycle, serial, scaleflags);
write_specgraph(cr, full, cycle, serial);
write_test_logs(cr, full, cycle, serial);
}
@@ -950,6 +955,7 @@ struct refine_args
int no_logs;
SymOpList *sym;
SymOpList *amb;
+ int scaleflags;
};
@@ -974,7 +980,7 @@ static void refine_image(void *task, int id)
do_pr_refine(cr, pargs->full, pargs->pmodel, pargs->verbose,
pargs->serial, pargs->cycle, write_logs,
- pargs->sym, pargs->amb);
+ pargs->sym, pargs->amb, pargs->scaleflags);
if ( crystal_get_user_flag(cr) == 0 ) {
pargs->prdata.refined = 1;
@@ -1013,7 +1019,7 @@ static void done_image(void *vqargs, void *task)
void refine_all(Crystal **crystals, int n_crystals,
RefList *full, int nthreads, PartialityModel pmodel,
int verbose, int cycle, int no_logs,
- SymOpList *sym, SymOpList *amb)
+ SymOpList *sym, SymOpList *amb, int scaleflags)
{
struct refine_args task_defaults;
struct queue_args qargs;
@@ -1027,6 +1033,7 @@ void refine_all(Crystal **crystals, int n_crystals,
task_defaults.no_logs = no_logs;
task_defaults.sym = sym;
task_defaults.amb = amb;
+ task_defaults.scaleflags = scaleflags;
qargs.task_defaults = task_defaults;
qargs.n_started = 0;
diff --git a/src/post-refinement.h b/src/post-refinement.h
index 8c345729..b8923d2c 100644
--- a/src/post-refinement.h
+++ b/src/post-refinement.h
@@ -61,10 +61,10 @@ extern const char *str_prflag(enum prflag flag);
extern void refine_all(Crystal **crystals, int n_crystals,
RefList *full, int nthreads, PartialityModel pmodel,
int verbose, int cycle, int no_logs,
- SymOpList *sym, SymOpList *amb);
+ SymOpList *sym, SymOpList *amb, int scaleflags);
extern void write_gridscan(Crystal *cr, const RefList *full,
- int cycle, int serial);
+ int cycle, int serial, int scaleflags);
extern void write_specgraph(Crystal *crystal, const RefList *full,
signed int cycle, int serial);
diff --git a/src/scaling.c b/src/scaling.c
index 39541926..5f896e65 100644
--- a/src/scaling.c
+++ b/src/scaling.c
@@ -389,7 +389,7 @@ static double total_log_r(Crystal **crystals, int n_crystals, RefList *full,
/* Perform iterative scaling, all the way to convergence */
-void scale_all(Crystal **crystals, int n_crystals, int nthreads)
+void scale_all(Crystal **crystals, int n_crystals, int nthreads, int no_Bscale)
{
struct scale_args task_defaults;
struct queue_args qargs;
@@ -448,10 +448,11 @@ void scale_all(Crystal **crystals, int n_crystals, int nthreads)
}
-/* Calculates G, by which list2 should be multiplied to fit list1 */
-int linear_scale(const RefList *list1, const RefList *list2, double *G,
- int complain_loudly)
+/* Calculates G and B, by which list2 should be multiplied to fit list1 */
+int scale_one(const RefList *list1, const RefList *list2, int flags,
+ double *G, double *B)
{
+ int complain_loudly = 0;
const Reflection *refl1;
const Reflection *refl2;
RefListIterator *iter;
@@ -474,6 +475,8 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G,
int n_part = 0;
int n_nom = 0;
+ *B = 0.0; /* FIXME */
+
x = malloc(max_n*sizeof(double));
w = malloc(max_n*sizeof(double));
y = malloc(max_n*sizeof(double));
@@ -556,6 +559,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G,
if ( r ) {
ERROR("Scaling failed.\n");
+ *G = 1.0;
return 1;
}
@@ -567,7 +571,7 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G,
int i;
for ( i=0; i<n; i++ ) {
STATUS("%i %e %e %e\n", i, x[i], y[i], w[n]);
- }
+ }
}
}
@@ -582,3 +586,16 @@ int linear_scale(const RefList *list1, const RefList *list2, double *G,
return 0;
}
+
+int scale_one_crystal(Crystal *cr, const RefList *list2, int flags)
+{
+ double G, B;
+ int r;
+
+ r = scale_one(crystal_get_reflections(cr), list2, flags, &G, &B);
+ if ( r ) return r;
+
+ crystal_set_osf(cr, G);
+ crystal_set_Bfac(cr, B);
+ return 0;
+}
diff --git a/src/scaling.h b/src/scaling.h
index 2161f69c..25aca31e 100644
--- a/src/scaling.h
+++ b/src/scaling.h
@@ -3,11 +3,11 @@
*
* Scaling
*
- * Copyright © 2012-2017 Deutsches Elektronen-Synchrotron DESY,
+ * Copyright © 2012-2018 Deutsches Elektronen-Synchrotron DESY,
* a research centre of the Helmholtz Association.
*
* Authors:
- * 2010-2017 Thomas White <taw@physics.org>
+ * 2010-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -38,12 +38,20 @@
#include "crystal.h"
#include "geometry.h"
+enum ScaleFlags
+{
+ SCALE_NO_B, /* Don't use Debye-Waller part */
+};
+
extern double log_residual(Crystal *cr, const RefList *full, int free,
int *pn_used, const char *filename);
-extern int linear_scale(const RefList *list1, const RefList *list2, double *G,
- int complain_loudly);
+extern int scale_one(const RefList *list1, const RefList *list2, int flags,
+ double *G, double *B);
+
+extern int scale_one_crystal(Crystal *cr, const RefList *list2, int flags);
-extern void scale_all(Crystal **crystals, int n_crystals, int nthreads);
+extern void scale_all(Crystal **crystals, int n_crystals, int nthreads,
+ int flags);
#endif /* SCALING_H */
diff --git a/tests/linear_scale_check.c b/tests/linear_scale_check.c
index 4cb6a1ba..5c723849 100644
--- a/tests/linear_scale_check.c
+++ b/tests/linear_scale_check.c
@@ -3,11 +3,11 @@
*
* Check that linear scaling works
*
- * Copyright © 2017 Deutsches Elektronen-Synchrotron DESY,
- * a research centre of the Helmholtz Association.
+ * Copyright © 2017-2018 Deutsches Elektronen-Synchrotron DESY,
+ * a research centre of the Helmholtz Association.
*
* Authors:
- * 2017 Thomas White <taw@physics.org>
+ * 2017-2018 Thomas White <taw@physics.org>
*
* This file is part of CrystFEL.
*
@@ -45,7 +45,7 @@ int main(int argc, char *argv[])
gsl_rng *rng;
RefList *list1;
RefList *list2;
- double G;
+ double G, B;
int r;
list1 = reflist_new();
@@ -70,7 +70,7 @@ int main(int argc, char *argv[])
set_partiality(refl2, 1.0);
}
- r = linear_scale(list1, list2, &G, 1);
+ r = scale_one(list1, list2, SCALE_NO_B, &G, &B);
STATUS("Scaling result: %i, G = %f\n", r, G);
return fail;
generated by cgit v1.2.3 (git 2.25.1) at 2024-11-16 04:43:45 +0000