aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--libcrystfel/src/geometry.h3
-rw-r--r--src/post-refinement.c115
2 files changed, 81 insertions, 37 deletions
diff --git a/libcrystfel/src/geometry.h b/libcrystfel/src/geometry.h
index 28d12033..b9fdf16d 100644
--- a/libcrystfel/src/geometry.h
+++ b/libcrystfel/src/geometry.h
@@ -81,7 +81,8 @@ enum gparam {
GPARAM_OSF, /* Linear scale factor */
GPARAM_BFAC, /* D-W scale factor */
GPARAM_ANG1, /* Out of plane rotation angles of crystal */
- GPARAM_ANG2
+ GPARAM_ANG2,
+ GPARAM_WAVELENGTH
};
diff --git a/src/post-refinement.c b/src/post-refinement.c
index bc449621..318d1733 100644
--- a/src/post-refinement.c
+++ b/src/post-refinement.c
@@ -200,38 +200,76 @@ static UnitCell *rotate_cell_xy(const UnitCell *cell, double ang1, double ang2)
}
+/* We set all the step sizes to 1, then scale them.
+ * This way, the size of the simplex stays meaningful and we possibly also
+ * avoid some roundoff errors */
+static double get_scale(enum gparam p)
+{
+ switch ( p ) {
+
+ case GPARAM_ANG1 : return deg2rad(0.01);
+ case GPARAM_ANG2 : return deg2rad(0.01);
+ case GPARAM_R : return 0.0005e9;
+ case GPARAM_WAVELENGTH : return 0.001e-10;
+
+ default : return 0.0;
+
+ }
+}
+
+
struct rf_priv {
const Crystal *cr;
const RefList *full;
enum gparam *rv;
int verbose;
+ const gsl_vector *initial;
};
-static void apply_parameters(const gsl_vector *v, enum gparam *rv, Crystal *cr)
+static double get_actual_val(const gsl_vector *v, const gsl_vector *initial,
+ enum gparam *rv, int i)
+{
+ return gsl_vector_get(v, i) * get_scale(rv[i])
+ + gsl_vector_get(initial, i);
+}
+
+
+static void apply_parameters(const gsl_vector *v, const gsl_vector *initial,
+ enum gparam *rv, Crystal *cr)
{
int i;
- double ang1, ang2, R;
+ double ang1, ang2, R, lambda;
UnitCell *cell;
/* Default parameters if not used in refinement */
ang1 = 0.0;
ang2 = 0.0;
R = crystal_get_profile_radius(cr);
+ lambda = crystal_get_image(cr)->lambda;
for ( i=0; i<v->size; i++ ) {
+
+ double val;
+
+ val = get_actual_val(v, initial, rv, i);
+
switch ( rv[i] ) {
case GPARAM_ANG1 :
- ang1 = gsl_vector_get(v, i);
+ ang1 = val;
break;
case GPARAM_ANG2 :
- ang2 = gsl_vector_get(v, i);
+ ang2 = val;
break;
case GPARAM_R :
- R = gsl_vector_get(v, i);
+ R = val;
+ break;
+
+ case GPARAM_WAVELENGTH :
+ lambda = val;
break;
default :
@@ -245,6 +283,7 @@ static void apply_parameters(const gsl_vector *v, enum gparam *rv, Crystal *cr)
crystal_set_cell(cr, cell);
crystal_set_profile_radius(cr, R);
+ crystal_get_image(cr)->lambda = lambda;
}
@@ -252,11 +291,14 @@ static double residual_f(const gsl_vector *v, void *pp)
{
struct rf_priv *pv = pp;
RefList *list;
+ struct image im;
Crystal *cr;
double res;
cr = crystal_copy(pv->cr);
- apply_parameters(v, pv->rv, cr);
+ im = *crystal_get_image(cr);
+ crystal_set_image(cr, &im);
+ apply_parameters(v, pv->initial, pv->rv, cr);
list = copy_reflist(crystal_get_reflections(cr));
crystal_set_reflections(cr, list);
@@ -288,6 +330,7 @@ static double get_initial_param(Crystal *cr, enum gparam p)
case GPARAM_ANG1 : return 0.0;
case GPARAM_ANG2 : return 0.0;
case GPARAM_R : return crystal_get_profile_radius(cr);
+ case GPARAM_WAVELENGTH : return crystal_get_image(cr)->lambda;
default: return 0.0;
@@ -295,26 +338,13 @@ static double get_initial_param(Crystal *cr, enum gparam p)
}
-static double get_stepsize(enum gparam p)
-{
- switch ( p ) {
-
- case GPARAM_ANG1 : return deg2rad(0.01);
- case GPARAM_ANG2 : return deg2rad(0.01);
- case GPARAM_R : return 0.00005e9;
-
- default : return 0.0;
-
- }
-}
-
-
static void do_pr_refine(Crystal *cr, const RefList *full,
PartialityModel pmodel, int verbose)
{
int i;
gsl_multimin_fminimizer *min;
- gsl_vector *v;
+ gsl_vector *initial;
+ gsl_vector *vals;
gsl_vector *step;
gsl_multimin_function f;
enum gparam rv[32];
@@ -333,7 +363,8 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
/* The parameters to be refined */
rv[n_params++] = GPARAM_ANG1;
rv[n_params++] = GPARAM_ANG2;
- //rv[n_params++] = GPARAM_R;
+ rv[n_params++] = GPARAM_R;
+ rv[n_params++] = GPARAM_WAVELENGTH;
residual_f_priv.cr = cr;
residual_f_priv.full = full;
@@ -343,17 +374,20 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
f.n = n_params;
f.params = &residual_f_priv;
- v = gsl_vector_alloc(n_params);
+ initial = gsl_vector_alloc(n_params);
+ vals = gsl_vector_alloc(n_params);
step = gsl_vector_alloc(n_params);
for ( i=0; i<n_params; i++ ) {
- gsl_vector_set(v, i, get_initial_param(cr, rv[i]));
- gsl_vector_set(step, i, get_stepsize(rv[i]));
+ gsl_vector_set(initial, i, get_initial_param(cr, rv[i]));
+ gsl_vector_set(vals, i, 0.0);
+ gsl_vector_set(step, i, 1.0);
}
+ residual_f_priv.initial = initial;
min = gsl_multimin_fminimizer_alloc(gsl_multimin_fminimizer_nmsimplex2,
n_params);
- gsl_multimin_fminimizer_set(min, &f, v, step);
+ gsl_multimin_fminimizer_set(min, &f, vals, step);
do {
n_iter++;
@@ -363,14 +397,22 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
if ( verbose ) {
double res = residual_f(min->x, &residual_f_priv);
- STATUS("%f %f residual = %e\n",
- rad2deg(gsl_vector_get(min->x, 0)),
- rad2deg(gsl_vector_get(min->x, 1)), res);
+ double size = gsl_multimin_fminimizer_size(min);
+ STATUS("%f %f %f %f ----> %f %f %e %f residual = %e size %f\n",
+ gsl_vector_get(min->x, 0),
+ gsl_vector_get(min->x, 1),
+ gsl_vector_get(min->x, 2),
+ gsl_vector_get(min->x, 3),
+ rad2deg(get_actual_val(min->x, initial, rv, 0)),
+ rad2deg(get_actual_val(min->x, initial, rv, 1)),
+ get_actual_val(min->x, initial, rv, 2)/1e9,
+ get_actual_val(min->x, initial, rv, 3)*1e10,
+ res, size);
}
- status = gsl_multimin_test_size(min->size, 1.0e-5);
+ status = gsl_multimin_test_size(min->size, 0.1);
- } while ( status == GSL_CONTINUE && n_iter < 100 );
+ } while ( status == GSL_CONTINUE && n_iter < 1000 );
if ( verbose ) {
STATUS("Done with refinement after %i iter\n", n_iter);
@@ -378,21 +420,21 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
}
/* FIXME: Not the right way to get the angles */
- ang1 = gsl_vector_get(min->x, 0);
- ang2 = gsl_vector_get(min->x, 1);
+ ang1 = get_actual_val(min->x, initial, rv, 0);
+ ang2 = get_actual_val(min->x, initial, rv, 1);
if ( rad2deg(fabs(ang1)+fabs(ang2)) > 5.0 ) {
ERROR("More than 5 degrees total rotation!\n");
residual_f_priv.verbose = 1;
double res = residual_f(min->x, &residual_f_priv);
STATUS("residual after rotation = %e\n", res);
residual_f_priv.verbose = 2;
- res = residual_f(v, &residual_f_priv);
+ res = residual_f(initial, &residual_f_priv);
STATUS("residual before rotation = %e\n", res);
return;
}
/* Apply the final shifts */
- apply_parameters(min->x, rv, cr);
+ apply_parameters(min->x, initial, rv, cr);
update_predictions(cr);
calculate_partialities(cr, PMODEL_XSPHERE);
@@ -403,7 +445,8 @@ static void do_pr_refine(Crystal *cr, const RefList *full,
}
gsl_multimin_fminimizer_free(min);
- gsl_vector_free(v);
+ gsl_vector_free(initial);
+ gsl_vector_free(vals);
gsl_vector_free(step);
}
generated by cgit v1.2.3 (git 2.25.1) at 2024-06-30 18:03:09 +0000